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ABSTRACT 

A consistent and stable finite-difference approximation of the 
vorticity transport equations has been applied to a plane Poiseuille 
flow with spatially periodic disturbances of finite amplitude. 

By application of the discrete Fourier transform to the solution of 
® Eoisson equation significant improvements in the speed and accuracy 
of the calculations have been obtained. Numerical tests indicate that 
the usual iterative schemes are inadequate if a large number of mesh 
points are employed. 

The effect of mesh size has been determined by extensive calcu- 
lations based on a consistent second-order representation of the linear 
eigenvalue problem. The results indicate that a fine mesh is required 
to accurately represent the behavior of even the large-scale unstable 
motions. The stability and accuracy of the present formulation was 
demonstrated, in the linear range, by comparison with the established 
results of the linear theory. 

Calculations in the non-linear range, at a moderate value of the 
Reynolds number, indicate the existence of an equilibrium state. The 
calculated values of the mean spectral density of kinetic energy are 
consistent with the exponential law predicted by Kraichnan for two- 
dimensional turbulence. The program, developed in this investigation, 
provides a means, of demonstrated accuracy, for the extensive investi- 
gation of two-dimensional turbulence in a shear flow. 
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I. INTRODUCTION 



A. BACKGROUND AND OBJECTIVES 

The phenomenon of turbulent flow has been the subject of numerous 
theoretical and experimental investigations since the publication of 
the results of the definitive experiments of Reynolds (1883). Yet, 
despite these efforts, our knowledge and ability to predict quantities 
of technical importance in a turbulent shear flow remain severely 
limited . 

The fundamental difficulty encountered in theoretical investigations 
of turbulence arises from the non-linear character of the equations of 
motion. The operation of statistical averaging introduces an additional 
dependent variable, the Reynolds stress, to the equations for the mean 
flow. The dynamic equation for the Reynolds stress involves higher 
order correlations, and so on. Although the basic equations are 
determinate, the statistically averaged equations are not. This is 
the so-called closure problem of the statistical theory. Hypotheses 
for the closure of the averaged equations range from the mixing length 
approach of Prandtl (1925) to the rather formidable Direct -Interact ion 
approximation of Kraichnan (1964). 

The advent of the high-speed digital computer has excited Interest 
in a new approach to the problem of turbulence, the numerical inte- 
gration of the basic equations by means of a suitable finite-difference 
representation. The potential value of such an approach lies not in 
the determination of the mean quantities in a particular flow, but 
rather in the possibility of determining, subject to the uncertainties 
which beset any numerical experiment, detailed and comprehensive 
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information about the structure of the turbulent motion. This infor- 
mation would be useful for the development of suitable models for the 
effects of the turbulence, and for the testing of existing closure 
hypotheses . 

The numerical approach is not without serious difficulties. The 
representation of the wide range of scales of motion involved in a 
turbulent flow produces staggering demands on the capacity and speed 
of the computer. In fact* the requirements are far beyond what is 
presently available (Orsag 1969). However, such calculations, at 
moderate values of the Reynolds number, nay be possible within a 
decade . 

A second difficulty concerns the development and testing of a finite- 
difference model which is stable, accurate and, in terms of computational 
requirements, efficient. A rigorous proof of the stability and conver- 
gence of finite-difference approximations to non-linear partial dif- 
ferential equations does not seem possible. Certain necessary conditions 
can be obtained by analysis of the linearized equations, but sufficient 
conditions are generally not available. The usual recourse is to 
demonstrate that the numerical solution does give a knewn result in 
some special case. A number of distinct finite-difference repre- 
sentations have been developed and employed for fluid dynamics problems, 
each with different strong points and weaknesses. There does not seem 
to be any general agreement as to which approach might be the best for 
a given application. 

The present capabilities of digital computers do permit numerical 
integrations for three-dimensional shear flows if a fairly coarse mesh 
is employed. In a number of recent investigations attempts have been 
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made to determine the behavior of the mean flow by means of a finite- 
difference representation which includes a model for the subgrid-scale 
motion (Deaidorff 1970), (Galloway 1968).' These are considered in 
somewhat greater detail in Section C.2. There are three fundamental 
difficulties associated with such undertakings. The first is the 
present lack of a sufficiently accurate model for the subgrid-scale 
motion. Secondly, the effect of the mesh size on the model for the 
subgrid-scale motion, and on the mean flow is largely undetermined. 
Finally the requirements for the accurate representation of even large- 
scale motions in the high shear layers near the boundaries may be 
greater than the coarse mesh permits. Despite these difficulties, 
qualitatively satisfying results have been obtained. 

Another approach, which is within the capabilities of current 
digital computers, is the numerical integration of the two-dimensional 
equations. There are fundamental and significant differences between 
two- and three-dimensional unsteady motions, which are considered 
briefly in Section B. However, the equally fundamental non-linearity 
of the equations remains. Two-dimensional disturbances in a parallel 
shear flow will, to a limited degree, excite motions of different 
scale and create apparent stresses which deform the mean flow. The 
results of a numerical simulation of a two-dimensional turbulent motion 
might be useful, to some degree, for the development of a suitable 
model for the effects of the small scales of motion, and the testing 
of various hypotheses concerning the structure of the turbulence. 

Considerations affecting the finite -difference formulation and 
solution are qualitatively the same for the two-dimensional problem as 
for the ultimate three‘ dimens iona 1 one. Whatever difficulties beset 
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the former will be increased in the latter, due to the greatly 
increased number of mesh points . Schemes , which are unsuitable in the 
two-dimensional case, will almost certainly be unsuitable in the three- 
dimensional one. Schemes which are stable and accurate in two dimensions 
are likely to be stable and accurate in three, although the fairly well- 
known stability constraints will be more severe in the three-dimensional 
case . 

For the case of a channel flow with periodic boundary conditions 
upstream and downstream, the two- dimensional problem offers the possi- 
bility of testing the stability and accuracy of the linear features of 
the formulation by comparison with the established results of the Orr- 
Sommerfe Id theory . 

The present investigation involves just such a numerical integration 
for a two-dimensional plane Poiseuille f]ow with periodic disturbances. 
The objectives are: 

1. To develop a stable, accurate and efficient finite-difference 
formulation for the plane Poiseuille flow, and to test its 
accuracy in a long-term integration by comparison with 
established results from the linear theory. 

2. To apply the method for the determination of finite-amplitude 
effects on the stability of small disturbances. 

3. To determine, if possible, an equilibrium condition for the 
two-dimensional flow by long-term integration with a growing 
disturbance . 

B. DIFFERENCES BETWEEN TWO AND THREE-DIMENSIONAL UNSTEADY MOTION 

There are fundamental differences between the kinds of non-linear 
interactions that can take place in two- and three-dimensional unsteady 
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flows. The effects of these differences are often described by noting 
the absence of "vortex stretching" in a two-dimensional flow. For 
zero viscosity, the circulation about all closed curves which move 
with a fluid particle remains constant. In a two-dimensional flow, 
all such curves remain in the plane, and the well-known result is the 
invariance of the vorticity of a particle. A three-dimensional flow 
is capable of bending or twisting a closed curve which moves with a 
particle. The curl of the velocity vector can change its direction: 
the component of vorticity in a specified plane can change. 

Kraichnan (1967) has shown that as a result of the existence of 
local inviscid constants, a two-dimensional flow admits of two formal 
inertial ranges, the usual -5/3 range and an additional -3 range. The 
-3 range, in this case, involves a cascade of vorticity to higher wave 
numbers, and the -5/3 range a cascade of energy to lower wave numbers, 
larger scales of motion. 

C. PREVIOUS NUMERICAL WORK 
1. Preliminary Remarks 

The equations of motion for an incompressible flow can be 
written in the familiar way, in terms of the velocity vector and the 
pressure . 

+ (v-v) V = - VP + W 2 v (1.1) 

with the additional requirement that 

v-v = 0 (1.2) 

to satisfy the conservation of mass. The equation for the pressure 
is found by taking the divergence of equation (1.1). 

V [ (V-V)V] = -V 2 P (1-3) 
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Another formulation, which implicitly satisfies the conservation 
of mass, is in terms of the curl of the velocity vector, F, and a vector 
velocity potential, \jf, which may be regarded as a generalized stream 
function . 

= - [(Vxi).V] f + <T-V) (Vxj) + W 2 r (1.4) 

where 

r = v(v-i) - v 2 | (i.5) 

and 

V = Vxi (1,6) 

Without restricting the solution, ^ can be specified to be non- 
divergent (Aziz 1967). Then (1.5) becomes 

r = - v 2 ^ (i.7) 

In the case of a two-dimensional motion, this formulation becomes the 
ordinary stream func tion-vorticity expression. 

Finite-difference approximations to the equations of motion in 
terms of V and P generally involve a staggered mesh arrangement in 
which pressures are defined at one family of mesh points and velocities 
at the other. The treatment of the boundary conditions for the pres- 
sure at the walls requires the establishment of virtual mesh points 
outside the domain and somewhat arbitrary specifications for the values 
of velocity at those points. Conservation of mass is not satisfied 
exactly, but the solution of a Poisson equation (1.3) for the pressure 
is designed so that at each point the divergence of the velocity is 
given a rate of change that will tend to correct for the accumulated 
divergence at that point. 
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Finite-difference representations based on equations (1.4) and 
(1.7) conserve mass implicitly. It is generally possible to derive a 
consistent explicit treatment of the boundary values for the vorticity 
without defining virtual points outside of no-slip boundaries. For 
three-dimensional integrations based on these equations a corrective 
procedure, similar to the one employed for the pressure calculation 
above, is required to maintain the zero divergence of \jr . 

The essential mathematical features cf both approaches are the 
same. At each time step a forecast for new values of V or T is 
obtained from the finite -difference representation of (1.1) or (1.4). 
The solution of a Poisson equation for the pressure in the first case, 
or three such equations for the components of \jr in the second, complete 
a time step. It is the Poisson solution which generally takes up most 
of the computing time. 

2 . Summary of Selected Numerical Investigations 

In a very well-known work, Fromm (1963) has employed a finite - 
difference representation based on the vorticity transport equations to 
follow the development of the vortex street behind a rectangular 
obstacle in a two-dimensional flow. The agreement between his computed 
streak lines and photographs of similar real flows is striking. The 
finite-difference representation involved an energy-conserving repre- 
sentation of the non-linear terms, an explicit formulation of boundary 
values of the vorticity at the moving walls, and periodic boundary 
conditions upstream and downstream. The time differencing was by the 
Du Fort-Franlcel method. The domain of numerical stability was deter- 
mined by analysis of the linearized equations and extensive numerical 
experiment. The finite -difference mesh contained 25 x 57 points. The 
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Poisson equation for the stream function was solved by Gauss-Seidel 
iteration . 

Perhaps the best-known f inite -difference representation based 
on the momentum equations is the Marker and Cell method of Harlow and 
Welsh (1965). The scheme, which was designed to be suitable for the 
representation of a flow with a free surface, employs a staggered 
mesh in which the components of velocity and the pressure are defined 
at distinct families of points. The treatment of solid boundaries 
requires the establishment of virtual points outside the domain, and 
somewhat arbitrary relations between values at the virtual and interior 
points. The failure of the representation of the non-linear terms to 
conserve kinetic energy is an important defect for reasons which will 
be considered in Chapter IV, Reported calculations employed first- 
order forward time differencing. The Poisson equation for the pressure! 
was solved by Gauss-Seidel iteration. 

An investigation which is closely related to the present work 
is due to Dixon (1966). It is an application of finite-difference 
methods to the problem of spatially-growing disturbances in Poiseuille 
and plane Poiseuille flow. The finite-difference equations are based 
on a version of the vorticity transport equations with an explicit 
separation of the laminar components. The representation of the non- 
linear terms failed to conserve kinetic energy. In addition the 
particular first-order explicit formulation of the boundary values of 
vorticity was inconsistent with the otherwise second-order finite- 
difference representation. Alternating direction methods were employed 
for the time differencing and the Poisson equation was solved by 
Successive Over-relaxation. For the plane Poiseuille flow Dixon 
assumed symmetry about the midstream, and employed a mesh with external 
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streamwise dimension 37 times the radius. The mesh contained 10 points 
between the centerline and the wall, and 296 in the direction of the 
stream. Conditions at the entrance region were s pec if ied, and the 
evolution of the disturbance was followed downstream. The results 
were somewhat uncertain. No evidence of growing disturbances was 
found. At very low values of the Reynolds number the disturbance 
decayed. In all other cases the behavior of the disturbance was 
character ized by periodic .amplification and decay. 

Pearson (1964) has considered several finite difference 
approaches to the solution of a test problem. He asserted that if 
the representation is based on the vorticity transport equations , an 
iterative method for the time step, which includes a method of fore- 
casting new values of the vorticity at the boundaries, is required. 

This was based on the observation that the basic equations do not 
include a specification of boundary values of the vorticity. It 
turned out, however, that the calculations were unstable unless the 
forecast was "smoothed" by weighting it with a fraction of the old 
value, and that the required weighting factor was in excess of 80 
per cent. .The procedure did produce a more rapid convergence of the 
iterative solution to the implicit equations for the time step. How- 
ever, the results do not seem to constitute a convincing argument for 
the necessity of an iterative treatment of vorticity at no-slip 
boundaries . 

Galloway (1968) has employed a second-order finite -difference 
representation of the momentum equations in an application to the 
three-dimensional flow in a channel of square cross section. The 
approach employs Du Fort-Frankel time differencing and an energy-conserving 
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representation of the non-linear terms. The Poisson equation for 

the pressure was solved by an interesting iterative procedure which 

takes advantage of the weak coupling between various groups of the 

algebraic equations. Calculations were performed in meshes containing 
3 3 3 

9 , 11 and 17 points. It was noted that the resolution was inade- 
quate to accurately represent the unsteady motion, however refinements 
of the mesh produced marked improvements in the results. An integral 
representation of the sub-grid-scale motion was proposed for subsequent 
calculations . 

Very recently Deardorff (1970) has performed calculations for 
a similar channel flow, at a high Reynolds number, which employed such 
an integral representation of the sub-grid-scale motion. The results 
are qualitatively satisfying, although it was noted that a number of 
improvements in the sub -grid model are needed. 
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II. STATEMENT OF THE PROBLEM AND BASIC EQUATIONS 



The equations of motion for a two-dimensional incompressible flow 



are 



5u Su chi §P 1 , c^u'N 

T _ +UT _ +V ^. .,_+_(_+ 

ox ay 



St 

Sv 



Sx 



Sy 

Sv 



Sx 



. Sv _ 

St U Sx V Sy Sy ^ R V-v 2 



Sp , 1 ^S 2 V , a 2 v'\ 

"N * T-k r\ • a ) 






(2.1) 

( 2 . 2 ) 



Sx Sy 

where u and v are the dimensionless components of velocity, P is the 
pressure and R is the Reynolds number which is defined by 



R = 



VL 



( 2 . 3 ) 



where V is the unit of velocity, L is the unit of length and v is 
the kinematic viscosity. The equation of continuity is 

s + < 2 -*> 

We consider the two-dimensional flow between fixed parallel planes, 
The origin of coordinates is midway between the planes with the x~axis 
parallel and the y-axis normal to them. 



t 



Let the reference unit of length be one-half the distance between 
the parallel planes, and the reference unit of velocity be the mean 
value of the steady laminar flow velocity distribution. 

Then the laminar flow velocity distribution and pressure gradient 
are given by 
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(2.5) 



U(y) = - f (y 2 - 1) 



dP = 3 * 

dx R 

If for the general case we let 



( 2 . 6 ) 









u (x,y ,t) 


= U(y) + 


u' (x,y, 


t) 






(2.7) 








v(x,y,t) 


= 


v ' (x,y, 


t) 






(2.8) 








p(xyt) 


= P(x) + 


p' (xyt) 








(2.9) 


Then equations 


(2. 


•1), (2.2) 


and (2.3) 


bee ome 










du' 

St 




‘ + 


, dU . , 

dx 


Su' , , 

&r + v 


Su' 

Sy ’ = ‘ 


Sp 

Sx 


Sp' , 

■ S7" + 


1 /'S 2 u' 

5 V 


*5 



(2-10) 



Sv ' , T , Sv ' 

sr + u sr 



3“' +|il =. o 



+ u * 



X, I 



dv ’ , , bv * = b? dp 1 1 f'cTv 1 , b^v 



dx + V dy 



dy dy R 2 
y ox 



2 "’ + <&--) 



V 

( 2 . 11 ) 

( 2 . 12 ) 



dx by 

The no-slip boundary conditions require that at y 3 + 1 

u 1 = v 1 - 0 

An equation for the disturbance pressure can be obtained by differenti- 
ating equation (2.10) with respect to x, equation (2.11) with respect 
to y, and summing. Thus 

' 2 ■ V*. 



5 p 1 , 5 p' 9 f du f dv 1 f dU L du bv H 

N 2 N 2 i _dx by vdy by J bx J 



bx by 



(2.13) 



We can, without confusion, delete the primes in the subsequent 
notation for the disturbance quantities. 

If we specify in addition that the disturbance makes no contribution 
to the total mass flow; 
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i.e. that 



(2.14) 



1 

f u(x,y,t) dy = 0 

-1 

then the reference unit of velocity, as previously defined, is as well 
the mean transport velocity. This conforms to the usual choice of 
reference velocity for turbulent channel flows. 

A more compact, and for purposes of computing, more efficient 
formulation of the basic equations is obtained by introducing the 
stream function and the vorticity. 



Let 



bjj Bty 
U ~ " By’ v ~ Bx 



(2.15) 



r - — — 

Bx By 



(2.16) 



Then equation (2.12) is implicitly satisfied. Equations (2.10) 
and (2.11) become 



where 



SC + „ £ . M.aiar + ai£ = irsfr + a?T) 

Bt Bx Bx By Bx Bx By R v^2 By^ 



r « ifi + ifi 

ay 2 



(2.17) 



(2.18) 



and 



U" = ~ = - 3 

dy 

The no-slip boundary conditions require that at y = + 1 



#-^=0 (2.19) 

ox dy 

Equation (2.14) becomes 
1 

/|dy = - Kx,-l,t) = 0 (2.20) 

-1 
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Thus the no-slip boundary conditions become 



i|r(x,+ l,t) = i|r(x,-l,t) = const 



( 2 . 21 ) 



(x 5 ±l,t) = 0. 



dy 



( 2 . 22 ) 



If we consider a finite region of the space between the parallel 



planes, and specify the initial state of the disturbance and the con- 
ditions at the upstream and downstream boundaries of the region, then, 
in principle, the evolution of the disturbance can be determined by 
means of a suitable numerical method. Unfortunately the question of 
the proper boundary conditions upstream and downstream is somewhat 
moot. A reasonable recourse is to assume that the motion is periodic 
in the direction of the mean flow. Thus we consider a disturbance 
such that 



where the basic period is 2h . 

It should be emphasized that, in the present formulation, the dis- 
turbance quantities do not describe a two-dimensional turbulent motion, 
but rather the total departure of the motion from the steady laminar 
flow. The turbulent motion is described by the difference between the 
disturbance and an appropriately defined mean value of the disturbance. 

In the present case, as a consequence of the assumed periodicity 
of the motion, the disturbance is homogeneous in the streamwise direc- 
tion. That is, all statistical correlations are independent of the 
origin of the x coordinate. Than an appropriate definition of the 
statistical mean value of a disturbance quantity is the average value 
over a period in the x direction. 



\|f(x + 2nh,y,t) = t(x,y,t) 
F (x ± 2 nh ,y , t ) = r(x,y,t) 



(2.23) 



(2.24) 
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Some remarks about notation are in order. It has become necessary 
to distinguish between quantities associated with the laminar flow, the 
disturbance, the mean flow and the turbulence. Some indulgence on the 
part of the reader will be required. 

With the exception of certain quadradic quantities such as turbu- 
lent energy or Reynolds stress, mean values, obtained by means of the 
linear operator 



are denoted by an overbar. All quantities associated with the two- 
dimensional turbulence are denoted by a circumflex. 

The distinction between local and mean values of quadradic quantities 
associated with the turbulence is indicated in the arguments of the vari- 
ous functions, or in the case of finite-difference expressions, the sub- 
scripts. Thus, since the averaging operator is linear, 



h 




dx 



-h 



U (x,y,t) = u (x,y,t) - u (y,t) 
v (x ,y , t ) = v (x,y,t) - v (y,t) 
(x,y,t) = (x , y , t ) - T (y,t) 

r (x , y , t ) = r (x , y , t ) - r (y,t) 



( 2 . 25 ) 



where 




( 2 . 26 ) 
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The mean flow velocity and vorticity are given by 



W (y,t) = u (y) - ^ I (y ,t) (2.27) 

Q (y,t) = - U" + f (y , t ) (2.28) 



It is convenient at this point to define a number of quantities 
which are calculated in the numerical scheme, and to derive expressions 
for them in terms of the stream function. The kinetic energy in the 
mean flow is 



E (y,t) = 1/2 w 2 (y,t) 

= 1/2 [u(y) - ^ ? (y,t)] 



(2.29) 



The normalized sum of mean flow kinetic energy over a domain which 
includes one period of the disturbance is 



h 1 

E(t) = 2h / dx \ f E ( y>t ) <ly 
-h -1 

1 

= \ j E(y>t) dy 

-1 



(2.30) 



The kinetic energy in the turbulent motion is 



* - s [(- 1/ + Ci)'. 



(2.31) 



The mean value is 

h 

E (y,t) = ~ J E (x,y,t) dx 
-h 



(2.32) 



And the normalized sum of turbulent kinetic energy over the domain is 
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( 2 . 33 ) 



h 1 

*E (t) = J dx ~ j E (x,y,t) dy 
-h -1 

1 

= \ f E (y,t) dy 

-1 

Finally the Reynolds stress produced by the two-dimensional turbulent 
motion is 



s (y,t) 



h 

1_ f M 

2h J ox a y 

-h 



(x,y,t) dx 



27 



III. ANALYSIS OF THE TRANSFORMED EQUATIONS 



A. BASIC EQUATIONS IN TERMS OF FOURIER COMPONENTS 

The equations of motion for the plane Poiseuille flow with periodic 
disturbances can also be expressed in terms of the wave number in the x 
direction by means of the finite Fourier Transform. 

h 

Let f n (y,t) = ^ elCXnX dx 

-h 

h < 3 - L) 

g n (y>t) = ^ f r (x,y,t) e ianx dx 
-h 



Then 



^ (*,y,t) 



s 



n=-co 



f n (y.t) e 



-ianx 



r (x,y,t) 



2 



n r=-co 



8n< y .O e 



-ianx 



(3.2) 



where a = n/h . Since both and F are real, 

f n (y>t) = f * (y»t) 

-n n 

g_ n (y>t) = s* ^ y,t ^ 



(3.3) 



where the asterisk indicates a complex conjugate. 
By means of the transform operator 

h 

1 f ianx , 

2h ) G dx 

-h 

and the associated orthogonality relation 
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h 

1 f -ia(n-s)x . . . 

2 ^ l e dx = 1 if s = n 

-h 

= 0 otherwise 

it can be shown that equations (2.17) and (2.18) are satisfied 

K - ° P ' a2 S„> - R S 8„ + ianR [(U-O) g n - ) f n ] 

- iaR H =0 
n 



if 



(3.4) 



n 



0 , 1 , 2 , 



where 




2 2 

a n 



f 

n 



(3.5) 



H = 



= £ { (n- 
s=l l 



s) 






n-s 



‘ 8 



f n-s> 



H- (n-rs) 



(fr 



8 n+s" 8 s' f n+s >J 



(3.6) 



and primes denote derivatives with respect y. The no-slip boundary 
conditions require that at y B + 1 

f = f = 0 (3.7) 

n n v 7 

n = 0 5 1 , 

The functions f (y,t) and g o () 7 ,t) represent the mean values of 
the disturbance stream function and vorticity respectively. In terms 
of these variables, the mean velocity, equation (2.27), is given by 



w(y ,t) = U (y) - f^(y,t) (3.8) 

The non-linear interaction between the n t ^ mode of the disturbance with 
the distorted mean flow is represented by the term 

(U-f *) g n - (U"-f"' ) f 
on on 

t h 

in equation (3.4). The non-linear interaction between the n mode 
and the other disturbance modes is given by H^. The equation for the 
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mean flow is, from equations (3.4), (3.5) and (3.6), 



£111 



o 



R Tt f 0 - l0R 1,1 s IX £ s> ■ ( S s f H ] " 0 < 3 ' 9 > 



s — 1 



Equations (3.4) and (3.7) govern the behavior of periodic disturb- 
ances of finite amplitude in a plane Roiseuille flow. They are an 
infinite set: of coupled non-linear partial differential equations. 

Previous analysis of the finite-amplitude disturbances has gen- 
erally involved the truncation of this set of equations, the deletion 
of certain non-linear terms to uncouple the remaining equations and 
various assumptions about the time dependence of the remaining modes. 

A brief description of the most prominent investigations is given in 
Section C. 

B. LINEAR THEORY 

If the amplitude of the disturbance is sufficiently small so that 
the terms quadradic in the amplitude, the non-linear terms, may be 
ignored in comparison with the linear terms, then equations (3.4) and 
(3.5) become 



In this case the solutions for the various modes are uncoupled. By 
absorbing the index n into the definition of CC, and assuming the sepa- 
ration of variables 




(3.10) 




n 



(3.11) 



f(y,t) = $(y) e ipt 



(3.12) 



g(y,t) = 6(y) e lpt 
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where 3 is, in general, complex. 



Equations (3.10) and (3.11) become 




(3.14) 



(3.13) 



Or in terms of $ alone, 




(3.15) 



The boundary conditions, equation (3.7), require that at y = + 1 



If the spatial wave number, a, and the Reynolds number, R, are 
specified, then equation (3.15) has solutions which satisfy the homo- 
geneous boundary conditions only for discrete values of the complex 
number 3* Equation (3.15) is the well-known Orr-S ommerfe Id equation. 

From equation (3.15) it is evident that, if the imaginary part of 
3 is negative, the amplitude of the disturbance increases with time. 

The disturbance is said to be unstable. If the imaginary part is 
positive the disturbance is stable. Of particular interest are con- 
ditions of neutral stability. The largest value of the Reynolds number 
at which, for all a, all the eigenvalues of equation (3.15) are stable 
is termed the critical Reynolds number. Schensted (I960) has shown 
that for a plane Poiseuille flow the eigenfunctions of equation (3.15) 
form a complete set. Any continuous function which satisfies the 
boundary conditions, equations (3.16), can be represented as a linear 
combination of them. Thus if, at given values of a and R, all of the 
eigensolut ions of equation (3.15) are stable, the flow is stable with 
respect to all infinitesmal disturbances . 



0 



(3.16) 
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The solutions of equation (3.15) for the case of plane Poiseuille 
flow have been extensively investigated. The zone of linear instability 
in the Reynolds number-wave number domain, and the associated curve of 
neutral stability, have been established in a number of invest igat ions , 
the most prominent of which is due to Lin (1945). Following a procedure 
based on asymptotic expansions, which was outlined by Heisenberg (1924), 
Lin calculated values of a and R on the neutral curve and determined a 
critical Reynolds number (in terms of the reference length and velocity 
of the present investigation) of 3540 at a wave number of 1.04. 

Thomas (1953) employed a basically fourth-order finite-difference 
representation of equation (3.15) for the determination of the primary 
eigenvalues associated with the symmetric modes at eighteen points in 
the vicinity of the critical point. The discrete representation of the 
eigenfunction involved 101 points between the centerline and the wall. 

The resulting algebraic equations were solved by direct Gauss elimination. 
The calculations, which were performed on an early electronic calculator, 
required two weeks of machine time. He determined by interpolation a 
critical value of the Reynolds number equal to 3850 at a wave number of 
1.04. For. the calculation at R = 6666.7 and a - 1.0 the values of the 
primary eigenfunction at the 101 grid points are tabulated. 

Potter (1966) has used methods generally similar to those of Lin, 
but with a different expansion for the inviscid solution, to investi- 
gate the linear modes of a combined plane Couet te-Poiseuille flow. 

His results for the Poiseuille flow indicate that there are slight 
errors in Lin's neutral curve for a > .8. Potter suggested that this 
was due to the slow convergence of Lin's expansion in that region. 
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More recently Grosch (1968) has employed an expansion of the 
solutions of equation (3.15) in terms of the known eigenfunctions 
of 

- a 2 $" + a 4 § = \ 4 § 

3 > (+ 1 ) = $' (± 1 ) = 0 

The expansion was truncated generally after thirty terms, and the 
resulting algebraic eigenvalue problem was solved by means of the 
Q R algorithm (Wilkinson 1965). His results showed excellent agree- 
ment with those of Thomas for both the various eigenvalues and the 
tabulated eigenvector at R = 6666.7, a = 1. The results for the 
neutral curve confirm Potter's observations. Solutions are obtained 
for both symmetric and antisymmetric modes. Instabilities are found 
only for the symmetric modes, and in no case is more than one of the 
calculated modes, at a given a and R, unstable. 

A portion of the neutral curve for the Plane Poiseuille flow based 
on the data of Lin and Grosch is shown in Fig. 1. 

All reference to Reynolds numbers and dimensionless wave numbers 
has been in terms of the reference length and velocity employed in 
this investigation, one half the channel width and the mean laminar 
velocity. In the investigations of Lin and Thomas the reference unit 
of velocity is the maximum laminar velocity. Based on that value, 
the Reynolds number is exactly 1.5 times values given here. In the 
works of Grosch and Potter (for the case of zero Couette component) 
the reference unit of length is the channel width, and the reference 
velocity is the maximum laminar velocity. Wave numbers are twice 
those given here, and Reynolds numbers three times those given here. 
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FIGURE I CURVE OF NEUTRAL STABILITY FOR PLANE PCSSEUiLLE FLOW. 



It is convenient, and useful for the interpretation of some of the 
results of the numerical integrations of the linearized equations, to 
indicate at this point certain relations for the total kinetic energy 
in the linear solutions. Equations (3.2) and (3.12) for the linear 
solutions are equivalent to 

'l f (x,y,t) = 2e $ R (y) cos (ax-S R t) + ^(y) sin (ax-g R t) (3.17) 



where 

3 = 3 R + i3j 

$ = $ & i$ 

k 1 



It is easy to show that, if h = vf a, the total disturbance energy, 
equation (2.33) is given by 



E(t) = a/2 






/ 

-1 



($■ 



R 



$j*) 



+ a 



<4 + 



4 i>] 



dy 



(3-18) 



_ -23 T t 

= const x e I 

The disturbance energy grows, or decays exponentially at twice the rate 
of the disturbance amplitude. The disturbance itself grows or decays 
exponentially and changes its phase at the rate |3 /a. 

If, in* the linear problem, the disturbance is composed of a number 
of the linear eigenfunctions with the same wave number, a, then the 
total disturbance energy will contain terms like equation (3.18) 
peculiar to each eigenfunction $ and eigenvalue (3, plus cross -coupling 
terms which involve pairs of the modes. 

For example, if the disturbance is composed of two eigenfunctions 
K x >y>t) = ^ (x ,y j t) 



where 



$(y) and 3 describe 



T](y) and 7 describe 



Then it can be shown that the total disturbance energy is given by 



E (t) = E L (t) + E 2 (t) 



+ CC e ^I^P^A cos [0 R "7 R ) fc j - B sin ^ (3 R "7 R ) t _J • 



(3.19) 



where E^(t) and E£ (t ) are from equation (3.18), and A and B are given 



by 



= / _ ($ R + *1 + ^ ($ R V + $ I J dy 



B “ / [ *R li • Ir) • ^ (, 1r *1 - *R ’ll)] d >’' 



-l 



Thus, in this case, the total disturbance energy contains terms 

with exponential growth rates - 2 p^, - 2 7 ^. and - ( 3 ^ + 7 ^), and the 

last term oscillates with frequency (3 - 7 ) , 

K R 

If 

- 23 I > -2 7][ 

-a 

then eventually E^(t) will overshadow the other terms.* 

C. NON-LINEAR ANALYSIS 

A fundamental consequence of the non-linearity of the equations of 
motion is the production of the Reynolds stress and the resulting de- 
formation of the mean flow. The distortion of the mean flow, which is 
dependent on the amplitude of the disturbance, can alter the rate of 
transfer of energy from the mean flow to the disturbance and, as a 
result, change its rate of growth or decay. 
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A second non-linear effect is the excitation of disturbance motions 



of other wave lengths. Dynamical constraints on the non-linear inter- 
actions in a two-dimensional motion limit the systematic transfer of 
kinetic energy to smaller scales of motion (Lilly 1965 , Kraichnan 1967). 
As a result, it is expected that, in terms of the kinetic energy in 
the various disturbance modes, the primary mode will remain dominant. 
Thus the question of non-linear stability centers on the behavior of 
the amplitude of the primary mode. If, with increasing time, the 
amplitude approaches zero, it is stable; if not, it is unstable. 

Stuart (1958) has employed the terms subcritical and supercritical to 
describe the stable and unstable regions in the wave number-Reynolds 
number domain for infinitesmal (linear) disturbances. Subcritical 
instability refers to the growth of a disturbance of finite amplitude 
at a wave number and Reynolds number for which it is stable in the 
linear analysis. 

A question separate from that of stability concerns the existence 
of an equilibrium state. It seems reasonable to expect that at some 
sufficiently large amplitude the non-linear interactions will limit the 
growth of the disturbance motion. 

In an early investigation, Meksyn and Stuart (1958) employed a 
variation of the asymptotic methods of the linear problem which 
included the distortion of the mean velocity. Their results suggest 
the existence of subcritical instabilities at values of the Reynolds 
number as low as 1933.3. 

Subsequently Stuart (1958) derived an approximate expression for 
the time dependence of the amplitude of the primary mode of the dis- 
turbance. The analysis is based on an energy balance and, as Stuart 
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indicates, may have questionable applicability to flows, such as the 
plane Poiseitille flow, for which the precise details of the mean motion 
are of crucial importance to the question of stability. In the develop- 
ment only the interaction between the primary mode and the mean flow is 
considered; the generation of higher harmonics in the disturbance 
motion is ignored. In addition it is assumed that the time derivative 
of the mean velocity is always negligible and that the phase velocity 
and the shape of the primary mode are as given by the linear theory. 

With reference to equations (3.4) this is equivalent to 
f n ( y jt) “ 0 n " 2 , 3 , 

= A(t) $ 1 (y) e lf ^ 

P R << 1 

where A(t) is an amplitude factor and $(y) and p are from the linear 
analysis. From the energy balance it is deduced that A(t) satisfies 
an equation of the form 



d(A 2 ) 

dt 



(k 1 -k 2 ) a 2 




where is associated with the energy production terms, with the 
viscous dissipation, and is determined by the non-linear interactions. 
At sufficiently small amplitudes the latter term is negligible and the 
sign of (K^-K^) determines the stability. If < 0, then regardless 
of the sign of (K^-K^) there is the possibility of finite amplitude 
instability. 

In subsequent undertaking, Stuart (I960) and Watson (1960) derived 
a formal expansion for the solutions to the non-linear equations (3.4). 
Based on certain order-of -magnitude arguments, the Fourier expansion 
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(3.2) for the solution is truncated after n=2 , and enough of the non- 
linear terms are deleted so that the remaining equations are uncoupled. 
A separation of variables is assumed for the functions f^(y,t) and 
f 2 (y,t), in which they are expressed in terms of perturbations of the 
corresponding linear solutions, and a complex amplitude function, A(t), 
which includes the changes in the phase speed due to the non-linear 
interactions. This function is assumed to satisfy 

diAL = _ Pi |a| 2 + a] |A| 4 + •■• (3.20) 

The result < 0 is taken to indicate either the existence of super- 
critical equilibrium condition or finite -amplitude stability in the 
subcritical domain. If a^ > 0 then, in the subcritical case an 
unstable equilibrium is implied, and in the supercritical case this 
result is taken to indicate the absence of an equilibrium condition 
(at least at small amplitudes). The last conclusion must be qualified 
since if the disturbance continues to grow then eventually higher order 
terms in (3.20) will become significant: the assumptions upon which 

the formal expansion is based will no longer apply. 

Reynolds and Potter (1967) have developed a somewhat more straight- 
forward formulation of the Stuart -Watson expansion and performed the 
indicated calculations for the case of a combined plane Couette- 
Poiseuille flow. The results for Poiseuille flow indicate values of 
a^ < 0 only near the lower branch of the neutral curve. 

In similar but more extensive calculations, Pekeris and Shkoller 
(1967) found the zone of negative a^ in the wave number-Reynolds number 
domain. It includes the lower branch, but excludes the vicinity of the 
upper branch of the curve of neutral stability for the linear problem. 
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IV. THE FINITE DIFFERENCE EQUATIONS AND SOLUTION 



A. CHOICE OF THE METHOD 

The formulation of the finite-difference model for the equations 
of motion is determined by considerations of accuracy, numerical 
stability and computational efficiency. The latter consideration 
refers to some weighted measure of computing time and storage 
requirement . 

The fundamental equations of motion contain three dependent vari- 
ables : pres sure and two components of velocity. The vorticity trans- 
port equations contain two: stream function and vorticity. The 
essential mathematical features of each are the same: non-linear 

parabolic equations for the time derivatives; a linear elliptic 
equation for the pressure in the first case, and the stream function 
in the second. Thus, although there are important differences in the 
treatment of the boundary conditions, similar finite difference 
methods can be applied to each set of equations. The same kinds of 
stability and accuracy constraints on the f inite -d i f ference parameters, 
6x and 6t, will occur in each case. 

The stability and efficiency of the f inite -d if ference representation 
of such equations is principally determined by the method of time dif- 
ferencing. A number of direct and implicit methods are available. 

Consider equations of the form 

= D(f,g); A(f ) = B(g) 

where 

f - f (x ,y , t) g = g(x,y,t) 
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and A, B and D are differential operators in x and y. 
representation in t let 



f(x,y,t K ) = f K (x,y) 
g(x,y,t K ) = g (x,y) 



For the discrete 



where 

t K - K6t 

Both Forward, and Modified Euler implicit time differencing 
require the simultaneous storage of two levels of the variables 
involved in the time differencing. 



N+l _ N 

6t 

r N+l r N 
f - r 

6 1 



N N 

D (f g ) + 0 (6 t ) 



1 ^ , N+l N+l. , 1 /£ N N. , n 

2 D ( f 8 ) + 2 » ( f S ) + 0 



(6t 2 ) 



If, in the case of Modified Euler differencing, the implicit 
N+l 

equations for f are solved by Gauss-Seidel iteration then it can 
be shown that the total storage requirement is the same as for forward 
d if ferencing . 

Central time differencing requires' the simultaneous storage of 
three levels of the variables involved in the time differencing. 



^N+l N- 1 
t - f 

26 1 



N N 2 

= D (£g) + 0 (5t ) 



Thus, for an M x N f inite -d iffcrence grid, the basic storage 
requirements may be summarized as follows. 
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TABLE I 



Computer Storage Requirements 





Momentum 

Equations 


Vorticity 
Transport 
Equat ions 


Forward, or Modified Euler 
Time Differencing 


5 (M x N) 


3 (M x N) 


Centra 1 

Time Differencing 


7 (M x N) 


4 (M x N) 



Clearly, in terms of the storage requirement, a formulation based 
on the vorticity transport equations is superior. In addition there 
are advantages in computing time. The velocity components, when they 
are required, can be obtained from the stream function with fewer 
operations than are required to perform the additional evaluation of 
a time derivative. A determination of the pressure ’field requires 
the additional solution of a Poisson equation. If this is required 
then some of these advantages disappear. 

Of the three methods of time differencing considered, the Modified 
Euler method is, in principle, the most accurate and, in particular 
for non-linear problems such as this one, the least susceptible to 
numerical instability. However it is implicit: the solution of a 

set of simultaneous equations is required at each time step. The 
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most practical means of solving the equations is by Gauss Seidel 
Iteration or Successive Over-relaxation. Each step in the iteration 
will require a solution of the Poisson equation for the stream func- 
tion; this demands all the operations that are involved in the course 
of a complete time step by an explicit method. Thus if an explicit 
method, of similar truncation error, can be made stable by reducing 
the time step by a factor less than the average number of iterations 
for the implicit case, it. is to be preferred. 

Forward time differencing is subject to linear stability con- 

2 

straints on the ratios 6t/6y and 6t/6y , which are associated respect- 
ively with the diffusion and linear advection terms. Early experiments 
with this method indicated that even if these constraints were observed 
numerical instability occurred when the non-linear terms became suf- 
ficiently large. 

If the ordinary second order finite difference representation of 
the diffusion term is employed, then the central time differencing is 
unstable at all values of 6t. This instability is removed by the 
Du Fort-Frankel (1953) approximation to the diffusion term. There 
still remains a constraint on the ratio 6t/5y, which is associated 
with the advection terms. Experiments were conducted with both the 
Central and the Modified Euler methods. At values of 6t small enough 
so that the Central differencing is stable, the Modified Euler method 
required from six to eight iterations: that is between six and eight 

times the computational effort to produce essentially the same results. 
Increasing the magnitude of the time step produced a predictable 
instability with the Central differencing and an increased number of 
iterations for convergence with the Modified Euler method. At the 
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stability boundary for the Central differencing the ratio 6t/5y is of 
the order unity. Thus the Modified Euler method offers advantages 
only in a domain where the truncation errors associated with 6t may 
be sufficiently large so as to obscure some transients in the solution. 
Accordingly, a formulation based on the vorticity transport equations, 
which employs the Du Fort-Frankel Central time differencing was used 
in this undertaking. 

It should be noted that an important group of methods of time dif- 
ferencing has not been considered in the preceeding remarks. These 
are the Alternating Direction Implicit methods of Peaceman and Rachford 
(1955) and Douglas and Rachford (1956) and their variants. In each 
case the time increment is performed in two steps. In the first step, 
an intermediate value of the dependent variable is obtained by express- 
ing certain of the spatial derivatives implicitly as in the Modified 
Euler method, and the others explicitly. In the second step the 
arrangement is reversed. Each step requires the solution of a set of 
simultaneous equations. The essential feature of the method is that 
these equations, in matrix form, are t r i -d iagona 1 : they may be solved 

directly and without storage penalty by a part icular ly ,s imple case of 
Gauss Elimination (Richtmyer 1967). In applications to the equations 
for incompressible flow, the variations of the methods are due pri- 
marily to the particular treatment of the advection terms. 

In the present investigation, the formulation of the upstream and 
downstream boundary conditions as well as the conservative finite dif- 
ference representation of the advection terms introduce additional 
non-zero diagonals to the matrices which define the equations for each 
step. The matrices are no longer tr i -d iagona 1 : the essential advan- 

tages of the alternating direction methods do not obtain. 
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B. THE FINITE DIFFERENCE EQUATIONS 



The computational domain, coordinates and the arrangement of the 
finite difference mesh are shown schematically in Figures (2) and (3). 
The mesh spacings are denoted by 6x and 6y, the value of the time step 
by 6t. The mesh contains M x N points: M points in the x direction 

and N points in the y direction. The external dimensions are 2h and 
2 respectively. 

The particular treatment of the periodic boundary conditions 
requires the definition of virtual points just beyond the upstream 
and downstream boundaries . The boundaries lie midway between these 
points and the adjacent interior columns of points . The mesh is 
arranged so that points lie on the no-slip boundaries. Thus there 
are M intervals in the x direction and N- 1 in the y direction. 



Let 

r K,J ” 1 

*K,J ‘ * 

Uj - U< yj ) 

where 

x= (K-l) 6x-h + 6x/2 K = 1 ,2 , • • *M 

K 

Yj = (J-l) 6y - 1 J = 1 ,2 , * ' *N 

t 

t = T, 6t 
l , n 

n=l 

and 

6x = 2h/M 

6y = 2/ (N-l) 



(4.1) 



(4.2) 
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FIGURE 2 COORDINATES AND DOMAIN 



□ 
□ 
□ 

K = 0 I 2 M-l M M+l 

O INTERIOR POINTS 

□ — VIRTUAL POINTS FOR PERIODIC BOUNDARY CONDITIONS 
X — NO -SUP BOUNDARY POINTS 

FIGURE 3 FINITE DIFFERENCE MESH 
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In the calculations M is always even, and N odd. Thus the origin of 

N+l 

the coordinates is at the mesh point K = M/2, J = — - — . 

The finite difference approximation t.o the vorticity transport 
equation is 



~K,J = - U ~K,J + U" |^K,J + J* T + ~ t 
6t J 6x 6x K, J R D 

K. , J 

< * 2 , 1 * 2,1 

p^ — Q. — j 6 — iLj^ j 

K,J . 2 k>J _ 2 K,J 

6x 6y 



(4.3) 

(4.4) 



6 5 2 

where - f — , — tc , etc., are the usual second-order central difference 

OX t 2 

6x t 

representations, D T is the Du Fort-Frankel approximation for the 

K , J 

l 

diffusion term and J is the conservative representation due to 

K , J 

Arakawa (1966). It may be regarded as a central difference approxi- 
mation to tiie mean of three distinct, but equivalent, expressions for 
the non-linear terms. 



J(t,r) 



5^ ST 

5y ck ck c)y 




In terms of the central difference operators 



(4.5) 



i,j - 3 ( ja k,j + <j + <j> 



where 
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(4.6) 




Each of the finite-difference forms, JA, JB and JC is a second order 
representation of the non-linear terms. The mean value has a number 
of important and desirable characteristics which are discussed in 
Section C . 3 . 

In detail equations (4.3), (4.4) and (4.6). 



r l+l pl-1 
K,J ~ K, J 
26t 



U 



J 



r 




\lj - ilr 

V K+1,J V K- 1 , J l 

26x ' J K,J 




(4.7) 




2 



J 



6y 



r, 



k 5 j 




(4.8) 
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J l = — ' 
K, J 3 
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(4.9) 






The periodic boundary conditions are approximated for equation 
(4.7) by assigning values of F and at the virtual points as follows, 



. I . t 

^0,J V M,J 



r l = r l 

0,J M, J 



Wi,j “ h,j 



r l = r l 

M+l, J M,J 



J - 1,2, • • -N 



J = 1,2, • ' -N 



0.10) 






Equation (3.8) is solved for ijr by means of a discrete Fourier trans- 

K j J 

form. The periodicity of T and ijr in the direction of the mean flow is 
implicit in the solution of that equation. 

The no slip boundary conditions, equations (2.21) and (2.22), 
require that 

4.1 ‘ *K,N = CO " St K-1.2.---M 



and that, at the walls, the normal derivative of the stream function be 
zero. The first condition is satisfied most conveniently by assigning 
a zero value to the stream function at the wall. The second is employed 
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to derive a boundary condition for the vorticity in equation (4.7). 
Since, at the wall, all derivatives of the stream function with respect 
to x are zero, the vorticity there is given by 



r 

I = — i- w. 
w L 

dy 

By means of Taylor series expansions , the interior values of the 
stream function can be expressed in terms of the stream function and 

its derivatives at the wall. Based on these expansions, finite- 

difference expressions for which satisfy 



W - 0 
dy 

can be derived. For example, at the lower wall (J=l), we have 

, , , e c^ r , 6y 2 <5^ i' , 6y^ , 6y^ 5% . 

*2 = h + 6 * $ w + — ^2 w + — ^3 w + 24 — ■ • ‘ 



(4.11) 



dy 

2 ^2 



dy 



+ 3 “h + 26 y ^ w+ 4 ^r2 W+ 8^6"“-3 w+ l 6 24“M + 

dy dy dy 

where, for simplicity, only the index J has been retained in the 
notation. Applying (4.11), we have, from the first expansion 

v3, 



, 3jl u . f djt 



w s. 2 

oy 



3 w 



w + ' 



(4.12) 



6y“ ~ ^ dy~ 

Equation (4.12) is the first-order expression for the boundary 
values of the vorticity. It can also be obtained by applying 
equation (4.8) at the wall with the assumption of a virtual point 
outside the domain to satisfy (4.11). Although equation (4.8) is, 
in the interior region, of second order, the application at the 
boundary is not. A second order formula can be obtained by combining 
the two expansions. 
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w + '• • • 



(4.13) 



r 



8 *2 ' *3 




W 




By including additional interior points third and fourth order formulas 



may be derived. They are 



r 



108 ^ - 27 + 4 - 85 ^ 




w + • * * (4.14) 



w 




r « 



576 - 216 + 64 - 9 - 415 ^ 




w 




(4.15) 



An integration based on equation (4.7) is not self-starting: 
another method is required for the first step. In addition the sup- 
pression of a certain weak instability associated with the central 
differencing requires the periodic averaging over two consecutive time 
steps and the performance of a step by another method. These points 
are considered in greater detail in Section C.3. The Modified Euler 
scheme is employed for these purposes. The associated finite differ- 
ence equation is 





pl-fl 2 pH'1, pH’1 

k+l,j" K , J r K-1,J 




(4.16) 




^K+l,j" ^K-l, J . ,t 

26x + J K,J 




where J is as given by equation (4.9). 
h , J 
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Finite-difference expressions for the various quantities which describe 
the mean flow and the two-dimensional turbulent motion follow in a 
straightforward way from equations (2.25) through (2.34). 

In terms of the discrete representation the relation for the mean 
value of the disturbance stream function 



(y,t) = J t (x,y,t) dx 
-h 



becomes 



*J = 2h 4,J 6x 



1 ^ l 

= - s * 

M K=1 K ’ J 



(4.17) 



Thus, in the discrete representation, the averaging operator simply 
implies the arithmetic mean value over the index K, or the zero mode 
in a discrete Fourier representation with respect to x (See Section 
D) . The stream function which describes the turbulent motion is 



•M _ . 1 *r l 

^'k,j ^k,j ’ b 

The velocity and vorticity of the mean flow are given by 



(4.18) 



w, = U. 



ft T l 

*J + 1 - ♦j-l 
26y 



Q t = 



- U" + r \ 

J 



(4.19) 

(4.20) 



Equation (2.29), for the mean flow kinetic energy, becomes 



_ l 1L 

E J 2 L U J 



T l T l 2 

b+i - 



■j-n 



25y 



(4.21) 
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The normalized sum is, by the trapezoidal rule 






(4.22) 



N-l 



N-l 



= f* V E l = V E l 

2 H=2 J N 1 J=2 J 



Similarly, equation (2.31) for the kinetic energy in the turbulent 
motion becomes, 



E l = 

K, J 2/.. 



-M 'A 2 *4 2 

~'' k,j-h~ ’k, j- r r v K+i,j" v k-i,j~ ! 
26y J _ 26x J 



(4.23) 



The associated mean value is 



M 

F 4 = 2 V 
J M K=1 KjJ 



(4.24) 



And the normalized sum over the domain is 



M ] -M 

E 



(4.25) 



Equation (2.34) for the Reynolds stress becomes 



S, = 



/y t >v. I 

, M , - A 

1 ^ J V K+1,J V K-1,J ^ K , J+ 1 



s s 

K=1 



{ 



♦ 

2§x 



1 1 

K, J 



26y 






(4.26) 



C. CONVERGENCE OF THE NUMERICAL SOLUTION 

A rigorous proof of the convergence of the numerical solution to 
the solution of the non-linear partial differential equation is not 
possible. Even if the non-linear terms are ignored in the analysis, 
certain difficulties remain. However, reasonable assurance of con- 
vergence can be obtained by analysis of the linear equations, and by 
numerical experiments. Finally, the convergence and stability of the 
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scheme, in the linear range, can be demonstrated by calculations 



which employ as initial conditions the eigenfunctions of the Orr- 
Sommerfeld equation . 



linear initial value problems, a finite difference approximation of 
the partial differential equation which is both consistent and stable, 
is necessarily convergent. A finite difference approximation is con- 
sistent with the differential equation if the associated truncation 
error tends to zero, as the magnitudes of the time and space incre- 
ments tend co zero. It is stable if the growth of errors in the 
numerical solution is bounded. 

The consistency of the finite-difference equations is considered 
in Section C.l. An analysis of the stability of an analogous linear 
difference equation is given in Section C.2. Finally, in Section C.3., 
the stability of the complete non-linear system is considered. 

1 . Cons is tency 



Richtmyer (1967) has shown that, for certain properly posed, 



The truncation error in a finite difference approximation is 



obtained by means of Taylor Series expansions. The results for 



typical terms in equations (4.3) and (4.4) are 







+ • • • 
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truncation error in equation (4.3) is 

6t 2 ra 3 r 



6t 2 / 6x 2 + 6yS a 2 r 



6x 



R V 
2 



2 2 

6x 6y 



; ) 
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L a t 3 
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6 1 6 x + 6y a T 
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- u j ^3 + u 3 

J ax J ax 3 



6x 6y 

,2 
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.. £± + Bi <Il + 3 ........ __ 

3 ax vax ax ay 



By 



a_ ^ ar <£ 



ax* 



at £ TJ\ 

ax axa y y 



ar a 3 t . i a 4 n 
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(4.27) 



Sy ax 

&y 2 r sr 



+ — - 1 

2R Sx 4j 



+ ^_ ^ . 3 b_ / at £r_ _ ar a 2 tj\ + i 

6 LSx ^3 Sy V.Sy SxSy "Sy "SxSyy Sx By 3 2 



aV 

2R ax 4J 



Thus the truncation error contains terms proportional to 
2 2 . . 2~ 



6t 

R 



~~' 2* ~ 2") ’ 6t2, 6x2 and 6y2 
6x 6y 



The first terms arises from the Du For t -Franke 1 approximation. 
Compatibility requires that 6t -» 0 faster than 6x or 6y. The 
truncation error in equation (4.4) is 



S2A Bjv + Bjlr 

12 * 4 12 -v 4 

ox oy 



(4.28) 
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Similarly the truncation error associated with the second-order 



formulation of the vorticity boundary condition is 



6 -v 4 

dy 



(4.29) 



It should be noted that the approximation for the boundary values 
of the vorticity involves values of the stream function at time level 
l. This approximation is correctly centered with respect to the time 
differencing which involves values at 1+1 and l-l. 

The over-all order of the approximation is 



n r 6t 2 Ax 2 + 5y 

L R V , 2 . 2 
ox oy 



and the scheme is consistent if 
LIM 
6x 



2 2 2 
+ 6t + 6x + 6y 



6y 

6t 



1 



0 



2 2 2 

6t_ r &2L + 6yl> = n 

R V A 2 A 2 J 

6x 6y 



(4.30) 



2 . Stability Analysis of an Analogous Linear Equation 

The numerical solution of the finite difference equations is 
stable if the round-off errors due to the finite decimal representation 
of real numbers are bounded in their growth. In the case of linear 
equations, the errors introduced at a given step are propagated by the 
same system of finite difference equations. In addition, the linearity 
permits, in the analysis, the superposition of random errors introduced 
at each step. For the purpose of determining what constraints, if any, 
must obtain for stability, it is only necessary to consider the ampli- 
fication, by the equations, of errors introduced at the initial or an 
arbitrary step. 
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Two approaches are in wide use. The first involves the deter- 
mination of the spectral radius of the matrix which defines the time 
step (Smith 1965). The finite -difference equations are expressed in 
the form 

E, - AE,^ (4.31) 

Where E^ is a vector whose components are spatial values of 
the dependent variable, (or in this analysis, the errors), at time 
step l. If the matrix A, of dimension q, has q linearly independent 
eigenvectors, r) g , then an arbitrary initial error can be expressed in 
terms of them. 

*0 = S % % 

S“1 

It follows from equation (3.39) that 



E , = AtE o 



Thus if v are the eigenvalues of A then 
s ° 



E, = L c v 1 n 
.. s s 's 
s = l 



t 



(4.32) 



Thus is bounded if the largest of the modulii of is less than 
unity. In some fairly simple cases, which generally involve two- level 
finite-difference formulas such as Forward or Modified Euler differ- 
encing, the bounds can be determined in a straightforward way by 
application of theorems due to Gerschgorin and Brauer (Smith 1965). 

If multi-level formulas, such as central time differencing are 
employed then the expression of the finite difference equations in the 
form of equation (4.31) is not so easy. In addition the determination 
of the eigenvalues or the spectral radius of the resulting matrix can 
be a formidable task. 
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A second approach which generally offers fewer difficulties, 
but which, unlike the matrix approach, does not take into account the 
effect of the boundary conditions is due to von Neumann. It was intro- 
duced in a paper by O f Brien, Hyman and Kaplan (1951), and is employed 
here. The initial errors are expressed in terms of their Fourier com- 
ponents. If none of the modes grow with increasing time the scheme is 
stable. Because the equations considered in the analysis are linear 
the various modes do not interact, and we need only to consider the 
behavior of an arbitrary mode. 

We consider a finite difference equation for the initial errors 
which is linear, but analogous to equation (4.7). 



pl+1 v* f ^ 

K, J ’ K, J K+1 S J ~ K- 1 , J , K,J+1 " K, J- 1 

~^Tt ' u o — + v o — — 




E l - 
K+1,J Vk/ 



E. 
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K,jy 



6x 



+ E 



t 

K- 1 , J 



(4.33) 



Let 



E, 



K , J-H 



6t 

r x U 0 6x 



X = 
x 



26t 

RSx^ 



Then equation (4.33) becomes 



Trl+1. V 1 " 1 ^ 1- TT l 
K,J + K , [ L K-1,J~ | 
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0.34) 



(4.35) 
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Let 
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(4.36) 



Where A are determined from the specified initial errors and 
nm 



=. e 



anm 6 t 



"nm 



For stability it is required that 



K I < 1 

1 J nm 1 — 



lim 



Since for the linear equation ve need only to consider an arbitrary 
combination of spatial modes, the substitution 

, . 1 _ t i 7 K i 'v J 

L K,J = 5 e X e ' y 

in equation (4.35) will reveal the necessary constraints. After 

l~l i 7 K i 7 J 

cancelling the common factor £ £ x e y and rearranging, 

equation (4.35) becomes 



j (^ x cos 7 ^ + cos 7 ) + i (r v sin 7 v + r_ sin 7 



x 



x 



x 



.)]? 



1-X 



= 0 



(4.37) 



1+X 

Which has roots 



- A {<’ 



£ = TT; \(\ cos 7 + X cos 7 ) + i (r sin 7 + r sin 7 ) 

* 1 x x y y x x y y 



(4.38) 



(A. cos 7 cos 7 ) + i(r sin 7 +r sin 7 ) 

x xy y x xy y J 



2 2 
+ (1-X ) 



Fromm (1963) has performed the analysis for the case 6 x = 6 y, which in 
the present notation means 
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r = r 
x y 



X * X 

x y 



X = 2X 



Because of the obvious difficulty in determining bounds on the 
modulus of ^ i- n equation (4.38), Fromm performed the analysis in two 
steps. In the first step only the diffusion terms are treated: the 

transport terms are assumed to be zero. In the second step the pro- 
cedure is reversed. The following is an extention of Fromm's analysis 
for the case 6x 6y« 

Let 
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1! 
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Equation (4.38) becomes 



§ = + ifi + /(atfiQ) 2 + 1-X 2 J (4.39) 



For diffusion with no transport 0=0. 
Equation (4.39) becomes 



5 = w {“ ± f 2 + 1_x2 J 



(4.40) 



Note that 



- X < m < X 

If the radical in (4.40) is real then \^\ reaches a maximum for 



2 2 

a) = X where 



F = 4- n + 1 

b 1 nj. i\ Z 1 



(X+l) 

If the radical is imaginary then 



60 



5 -i+r{ m ± 1 ^ 2 "“ 2 ' 1 .] 



15! 



4 2 -i = jx ^ 7 

i+x h+i 



In each case |£| < 1 for all values of X and X . Thus, in 

x y 5 

relation to the diffusion terms, the Du Fort-Frankel Central differ- 
encing is stable at all values of 6t. This result is given by 
Richtinyer (1967) for the one-dimensional case. 

For transport with no diffusion, 0 ) = 0, and A. = 0. 

Equation (4.39) becomes 

§ = in + /- n 2 + i (4.4D 

2 

If Q <1 the radical is real and 

m = i 

2 

If, however, fi > 1, the radical is imaginary and 

§ = i n + Jq 2 - ij 

15 1 = 2n 2 - i + 2n 



For the second root 

m > i 

Thus, in relation to the transport terms the Central differencing is 
stable if 

2 2 
0 = (r sin 7 -f* r sin 7 ) <1 

x 7 x y y — 

Then from equations (4.34) 
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6t 



!» + !o 

6x 6y 



< 1 



A more conservative form is 



/ I u n I l v nk 

6t l — H — < 1 

v 6x by / ~ 



which is equivalent to 



6t < 



6x6y 



i u o I 6y + l v o |6x 



(4.42) 



Although a complete analysis of even the linear equations* is not 
carried out it is reasonable to expect that a constraint of the form of 
equation (4.42) will be required. 

3 . Stability and Accuracy of the Complete Non-linear Equations 
From the stability analysis of the linear equations it was 
expected that a constraint of the form of equation (4.42) would be 
required for the general case where both diffusion and advection terms 
are present, the latter including the non-linear terms. Accordingly, 
the value of the time increment was determined by 



6t 



f 6x 6y 

u o 6y + v o 6x 



(4.43) 



where u^ and v^ are the maximum absolute values of the x and y com- 
ponents of total velocity over a large sample of the mesh points, and 
f is a parameter which was determined by numerical experiment. The 
values of u^ and v^ were recomputed at regular intervals, and &t was 
adjusted accordingly. Experiments with a variety of mesh sizes and 
disturbance amplitudes indicated critical values of the parameter f in 
the range .95 < f^ < 1.02. When the critical value of f was exceeded, 
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the instabilities developed quite rapidly, always within 100 time 
steps and in most cases within 50 steps. At low disturbance ampli- 
tudes the instability developed first in the midstream region, where 
the laminar transport velocity is the greatest. At higher amplitudes 
of the disturbance, the instability occurred first in the wall region 
where, owing to the character of the disturbance, the advection was 
the greatest. The instabilities were, in most cases, unmistakable: 
an explosive growth in the numerical values of the disturbance 
quantities immediately preceded an abnormal termination of the calcu- 
lation due to numerical overflow. At values of f very close to the 
critical, a more insidious situation occurred in which the initial 
amplification of errors in the disturbance stream function and 
vorticity, and hence in the computed quantities u^ and v^, was in- 
hibited by the ensuing automatic reduction of 6t. The calculations 
were characterized by periodic amplification and decay of large ei'rors . 
The results were in poor agreement with recalculations with a reduced 
value of 6t. If the automatic recalculation of the time step was by- 
passed then the explosive instability described above would occur. 
Accordingly f was limited to values less than .6, where no such dif- 
ficulties were encountered. 

In similar numerical experiments with the Du Fort-Frankel 
Central differencing applied to the non-linear vorticity transport 
equations, Fromm determined a critical value of f equal to .8. An 
important difference between the finite-difference equations employed 
by Fromm and those in the present work is associated with the treat- 
ment of the non-linear terms. Fromm employed a representation of 
those terms equivalent to JB T of equations (3.6). As has been 

j J 
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noted, in the present case the mean of JA _ , JB , and JC is 

K y U XV y U XV y J 

employed. The significance of the various finite difference repre- 
sentations of the non-linear terms is considered below. 

Phillips (1959) has shown the existence of numerical instabilities, 
associated with the non-linear terms, which result from the discrete 
representation of continuous functions, but which, in principle, can- 
not always be eliminated by reducing the time step or the mesh spacing. 

This kind of instability is the result of uncontrolled "aliasing" or 
confusion of frequencies (Hamming 1962). To illustrate, note that at 
a given value of the index J the M values 



%J 



K = 1,2, * • -M 



M 

can be expressed in terms of ~ + 1 distinct Fourier modes. 

1 f . : 2nrK , , . 2rmK x , •) 

t T = «' a n T +2 ^ V a T cos — M 1- b sin — — j + a . T L 

Y K,J M 0,J V n,J M n , J M s M/2,J ( 

where the period associated with the highest mode number is 25x. 
Since , 



2nnK 2 iittK 1 2(n+m)TTK . 2(n-m)TiK 

cos — — cos — — = 77 cos r— z h cos ~ N ■ ■ 

M M 2 M M 



etc . , 



the non-linear interaction between modes n and m will produce excitation 
in modes n-hn and n-m. If nH-m is greater than M/2 that interaction will 

be misrepresented by the discrete system. 

M M 

Let m+n = 2 '+r, r < ~ 
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Thus the aliasing involves a spurious interaction with mode 
r = n-hn - M/2. If both m and n are less than M/4, that is if each 
wave has a period greater than 46x, no misrepresentation will occur. 

The instability that results from aliasing occurs after a large number 
of time steps and is characterized by explosive growth of the energy 
in the system. Phillips was able to suppress this instability by, 
"periodically eliminating all Fourier components with wave lengths 
smaller than 46x", that is, by zeroing the upper half of the repre- 
sentable spectrum. In effect, in each spatial or spectral coordinate, 
twice as many degrees of freedom are carried in the calculations as 
are retained in the solution. For calculations explicitly in terms of 
the spectral form of the equations this seems to be the only recourse. 

Arakawa (1966) has devised a number of finite difference forms 
for the non-linear advection terms which eliminate this instability 
by controlling or limiting the aliasing. They are various linear com- 
binations of the forms JA, JB and JC as given, in terms of the central 
difference operators, by equations (4.6). For a confined flow or one 
with identical conditions at inflow and outflow boundaries it is easy 
to show that, in the absence of viscosity, the spatial integrals of 
kinetic energy and any function of vorticity are fixed quantities. In 
a viscous flow the advection terms make no contribution to the rates 
of change of these quantities. Lilly (1965) has shown that the feature 
essential to the control of the aliasing is the conservation of mean 
kinetic energy or mean squared vorticity implicit in the finite- 
difference representation of the advection terms. A number of stable 
conservative schemes are given by both Arakawa and Lilly. The only 
scheme which conserves both mean kinetic energy and squared vorticity 
is 
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J(i|f,r) = J (JA + JB + JC) 

which is the form used in the present undertaking. 

There is, in addition, a third type of instability which is 
associated with the central time dif fereneing . This weak instability 
is a computational mode associated with the three-time level central 
difference representation of = f(F). In effect, the central dif- 
ference representation replaces a first order (in time) continuous 
equation with a second order discrete one, which has an additional, 
non-physical, solution (Lilly 1965). It takes the form of a slowly 
growing oscillation, of period 26 t, about the true solution. If it 
remains unchecked then eventually the solutions at even and odd time 
steps diverge. This computational mode can be suppressed quite simply 
by periodically averaging the solutions over successive time steps and 
restarting the calculations with the averaged values (Will iams ' 1969) . 

This periodic interruption of the sequence of the calculation 
provides a convenient opportunity for the readjustment of the value of 
6t. In the midst of the continued central time differencing 6t must 
remain fixed. However when the averaging takes place it can be 
adjusted in accordance with equation (4.43), a single time step per- 
formed by the Modified Euler method (which does not require values at 
the previous time level), and the normal sequence can be continued. 
This procedure proved to be quite effective for both purposes. No 
evidence of the growth of the computational mode was observed. 

Fischer (1965) tested the accuracy of four methods of time 
differencing on a one -dimens ional linearized version of the equations 
of motion, and found that the central differencing was decidedly 
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superior. In similar tests on the wave equation Kurihara (1965) 
found that, provided the computational mode is suppressed, it was 
superior in terms of measured errors in amplitude and phase velocity 
to the other methods tested. 



D. SOLUTION OF POISSON'S EQUATION FOR THE STREAM FUNCTION 

The finite-difference relation between the stream function and 
the vorticity is given by equation (4.8). In the course of each time 
step new values of the vorticity are obtained from equation (4.7). 
Equation (4.8) must then be solved for tt e matching values of the. 
stream function. 



r = ^k+UJ 

K,J 



- 2 Vj + 



iLdL-J + 



V J+i ~ 2 Vj 



+ *K. 



J-l 



(4.8) 



6x 



6y 



K = 1,2, • • -M 
J = 2,3, • • • N- 1 

In matrix form equations (4.8) are 

f = Bt (4.44) 

where F and ijr are vectors of dimension d = (N-2) M and B is a banded 
symmetric matrix of dimension d x d. The solution to equation (4.44), 

♦ - B _1 f 



is obtained in three steps by means of the discrete Fourier transform. 
In the transformed domain the finite-difference relation between the 
stream function and the vorticity has the form 



r 

n 



= A * 
n n 



n = 1,2 , * 




(4.45) 



where F and \[r are vectors of dimension N-2 and the A are symmetric 
n n n n 

tridiagonal matrices of dimension (N-2) x (N-2). Because the matrices 
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A are tridiagonal, the solutions to equations (4.45) can be obtained 

n 

by a particularly simple case of Gauss elimination (Richtmyer 1967). 
The discrete Fourier transform pair is given by 



T 

n 



M 

Z 

K=1 



F e 
K 



i 0 



n,K 



f k = 



M 



M 

X 

n-1 



T e 
n 



-i 9 



n,K 



n = 1,2 , • • -M 
K = 1,2, •• -M 



(4.46) 



Where T and are, in general, complex and, 



e = 2n(n-l) (K-l) 
n,K M 



<4.47) 



The associated orthogonality relation is 



M 

r- E e^ ^n ,s-K+l = 1 if s=k 

n=l 

= 0 otherwise 

For the solution of equation (4.8) let 



(4.48) 



M 



i 0 



G = £ T e n> K 

n > K=1 5 K = 1,2 , • • -M 



M 



i 0 



J = 1 ,2 , • • • N 



F = E T e n,K 
n , J ^ k , J 



(4.4 9a ) 



Then from equations (3.65) T v and \|r , are given by 

K , J K , J 



M 

1 -10 
r ir T = T 7 S G T e n,K 

K,J M n=l n,J K = 1,2, • ' -M 

, 1 “ -i9 K 1 - 2 >'-' N 

\ F n,j e n ’ K 
n=l 



(4.49b) 



From equations (3.1) and (3.2) it is evident that 



n,K 



(x R -c) 2n (n-1) 
2h 
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where 



c 



6x 

2 



h 



Therefore equations (4.48) and (4.49) implicitly satisfy the periodic 
boundary conditions on F and i|r . The no-slip boundary conditions 
require that 



F = F = 0 n = 1,2, • * *M 

n, 1 n,N ’ 5 

Substitution of equations (4.49) in equation (4.8) yields, 



M 

S 1 
n= 1 i 



6y 2 G - Ff tj-i" a F t + F tj .1 
n,J i.. n,J+l n n,J n,J+] 



e_1 9 n,K = 



0 



(4.50) 



K = 1,2, •• -M 
J = 2,3, • ' -N-l 

where 

“n " 2 L 1 •— 2' (COS e „,2 ‘'J < 4 - 51 > 



Equation (4.50) can be satisfied at all values of K only if 




C T = F - a F + F TJL1 

n , J n , J+ 1 n n , J n , J+ 1 



n = 1,2, • 


• -M 


J = 2,3,- 


• ‘N-l 



(4.52) 



Although F and G . are complex, the coefficients in equation (4.52) 
n , J n , j 

are all real. The solutions for the real and imaginary parts may be 

obtained separately. Because F , and \[? , are real, the real and 

K , J K , J 

imaginary parts of G , and F , are symmetric and anti -symmetric 

n , J n , J 

M 

respectively about n = ^ + 1. (Recall that M, as defined, is always 
even.) Thus M x (N-2) quantities in the x domain beget exactly 
M x (N-2) distinct quantities in the wave number domain. There are 
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exactly M distinct sets of equations of the form equation (4.52) in 
terms of real variables. 

Equations (4.8) can be solved in this manner without the heavy 
storage penalties that are usually associated with direct methods. In 
this case the additional storage requirements are trivial: an M 

dimensional complex vector for the row by row transformation of T 

K , J 

and the inverse transformation of \|f , and three N-2 dimensional 

k, J 

vectors to accommodate the solution of the tridiagonal systems. The 
vorticity array is transformed row by row, and the transformed 
quantities are stored in the array reserved for the stream function. 

The tridiagonal equations (4.52) are solved column by column. In each 
step, the solution for the transform of the stream function replaces 
the stored values of the transform of the vorticity which are no longer 
required. Finally the transform of the stream function is inverted 
row by row. In order to minimize the accumulation of roundoff errors 
in the direct solution, the calculations associated witli both the 
discrete transform and the Gauss elimination are performed in double 
precision. The equivalent of sixteen decimal digits are carried in 
the calculations . 

The transforms and inverse transforms are obtained by means of the 
remarkably efficient Fast Fourier Transform algorithm of Cooley and 
Tukey (1965). The computational effort required to perform the trans- 
form on an M dimensional vector is proportional to M log (M) . There 
are roughly 2 N such transforms required in the solution of equation 
(4.8). The effort required to solve each of the M tridiagonal systems 
is proportional to 3 N. Hence the total computational effort is pro- 
portional to NM (2 log (M) i 3). 
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In an early version of the program equation (4.8) was solved by 
Successive Over-relaxation. It is an iterative method which is easy 
to program, involves no storage pena It ies. and which, at least for 
small arrays, is rapidly convergent. 

Equation (4.8) can be written 



2 2 

V, J = , 2, 2*7 LW,J + +-1,J + ~2 0k,J+ 1 + + " 6x r K,J-i 

2 (ox +6y ) 6y 



(4.53) 



Equation (4.53) can be used as the basis for various iterative 
formulations. If the superscript denotes the iteration number then 



♦ 



i+i 



6v 



r , i , , t , 6x^ 

1 T + V . 



K,J 0/£ 2 2, L V K+1,J' y K-1,J‘ . 2 Ok.J+I' y K,J-U 

2 (6x +6y ) 5 6y 



i+ i'l 






6x2 r K.J 



(4.54) 

defines a Jacobi iteration. With this method separate computer storage 
must be reserved for the l and 1+1 iterat.es. If the current values are 
employed in the formula as soon as they are available then the extra 
storage requirement is eliminated. In addition the convergence rate 
is, for the Poisson equation, approximately doubled (Varga 1962). 
Equation (4.54) becomes 



♦ 



6y 



1+1 



K, J ... 2,. 2. 
2(6x +6y ) 



,t ,.l+l ,6x Al . , l+l ^ . 2 

. V K+1,J +V K-1,J + 6y 2 V V K,J+1 + V K,J-iy " 6x r K,J.J 



(4.55) 



Equation (4.55) defines a Gauss-Seidel iteration. 

By the addition of an acceleration or over-relaxation parameter to 
equation (4.55) the convergence rate can be considerably increased 



,1+1 ,, x ,1 . o_'6y r ,l , , 11*1 .ox [ji , iTi > 

+ ,J “ (1 "“ ) V K,J + 2(6x 2 + ^ y 2 ) L V K+1,J + Vl,/ fv 2 V V K,J+l +1; K,J-iy~ 6x r K,J-i 



t+1 4-^Vl/ +,!r t+1 \^x 2 l 

6/ 

(4.56) 
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Equation (4.56) defines the Successive Over-relaxation. In this 
case the convergence rate is strongly dependent on the choice of to. 

The optimum value of cd, which for the Poisson equation falls between 
one and two. is determined by the modulus of the largest eigenvalue 
of the matrix which defines the system of equations (Varga 1962), 
however it is generally easier to determine it by experiment. 

With each of the various iterative methods an approximate solution 
to the system of equations is obtained by s)rs temat ica 1 ly applying the 
iteration formula at each point in the field until some measure of the 
changes in the various values falls below a specified tolerance. In 
order to equalize the influence of the boundary conditions it is 
necessary to start successive sweeps at opposite boundaries. 

Early experiments with certain contrived test problems indicated 
that with the optimum value of CD, convergence was generally obtained 
with one fifth the number of iterations required for the Gauss -Seidel 
Iteration, and that the errors in the solution were reduced by as much 
as a factor of ten. Unfortunately in application to the present prob- 
lem the performance of even the Over-relaxation is inadequate. The 
computational effort required to perform one iteration in an M by N 
array is proportional to M N. The number of iterations required for 

1/4 

convergence to a fixed tolerance was roughly proportional to (M N) 

Finally the tolerance required for a fixed error in the solution varied 

1/2 

nearly inversely as (M N) . The over-all computational effort was 

7/4 

roughly proportional to (M N) . In the integration of the vorticity 
transport equations the relaxation solution to the Poisson equation 
required more than 80 per cent of the total computing time. More 
sophisticated relaxation techniques, such as line or block over-relaxation, 
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might double the convergence rate (Varga 1962), but the essential 
difficulties remain for sufficiently large dimensions. 

A number of numerical tests were conducted in order to compare 
the efficiency and accuracy of the various schemes. In these tests 

values of were established from 

i'w , J 



'K, J 






Y 2 + 

J 



6 

S 

m=l 



1/m 



s in 



(mx K ) 



v 2 

Y J C0S 



(mx K )_ 



Values of I* , were obtained from the prescribed values of by 

i) J K J J 

application of equation (4.8) , Thus the stored values of , are 

K , J 

an exact solution of the finite-difference equations (4.8) corres- 
ponding to the values of j. Solutions of equations (4.8) denoted 

by s , were obtained by Gauss-Seidel iteration, Successive Over- 
K , J 

relaxation and by the direct method. 

The mean and maximum relative errors defined by, 



E 

mean 



1 

M(N-2) 



N- 1 M 
E E 
J=2 K=1 



K, J 



1 K, J 



f K,J 



E 

max 



max 



l 



S K,J ~ ^K, J 



^K, J 



K = 
J = 



1.2, - • -M 

2.3, - • -N-l 



} 



were calculated for each test. 

The iterative processes were assumed to have converged when at 
every point in the field 



j l+l l 

l S K,J ~ S K,J „ 

s 

where s is the mean absolute value of s T , at the start of the 

K , J 

iteration procedure and g is an assigned tolerance. In order to 
minimize the computer time expended in testing for convergence, the 
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test was applied only at even-numbered sweeps. In addition the test 

was bypassed at all points following the first point at which it was 

not met. Initial values of s^ were established reasonably close 

to the exact solution in order to approximate the situation in the 

numerical integration, where values of at time t are used as a 

K , J 

first guess for values at t + 5t in the iterative procedure. These 
initial values were 



K, J 



+ Ti 



K, J 10 L' 



J 1 J L 



sin (2 x k > + cos (5 x k )J 



In every case the experimentally determined optimum value of the 
relaxation parameter, cr>, was employed for the Successive Over- 
re laxat ion . 

Internal timers were started and stopped at the beginning and at 
the end of each solution and the elapsed times, in seconds, were 
printed along with errors. The results are given in Table II. They 
indicate rather dramatically the superiority of the direct solution in 
terms of both accuracy and speed. 

Tests with the iterative methods were performed at two and three 
values of the convergence tolerance in order to verify that the pro- 
cedures were still converging systematically at the specified tolerances. 
Generally a reduction of the convergence tolerance by a factor of two 
produced a solution with the errors reduced by about the same factor. 

The effects of the rapid reduction of the convergence rates of the 
iterative methods due to increased array dimensions are quite evident. 

The Gauss-Seidel iteration is, in particular, quite useless in 

application to large arrays. For the 64 by 201 point mesh, a compari- 

— 4 — 5 

son of the tests at g = 10 and 5 x 10 shows that the expenditure of 
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Test Solutions of Poisson Equation 



eight minutes of computer time, in the performance of 522 iterations, 
reduced the mean error by only twenty-five percent, and that the 
errors in the resulting solution are such that it contains no signifi- 
cant figures! 

The performance of the Successive Over-relaxation is considerably 
better, though still inadequate for the long-term integration. The 
results for the 128 by 101 and the 64 by 201 point arrays indicate 
that solutions with two significant figures can be obtained in from 
two to foul' minutes of computer time. 

The rate of convergence which occurs at any point in an iterative 
solution is proportional to the difference between the current esti- 
mate and the solution. As the gap is reduced the rate of approach 
decreases to some limiting value. It is this asymptotic convergence 
rate which is the significant factor in the determination of the 
computational effort required to obtain a solution of specified 
accuracy (Varga 1962). The three tests with the Successive Over- 
relaxation in the 64 by 201 point array permit the estimation of two 
successive convergence rates. The second rate is only very slightly 
less than the first: their ratio is about .97. These are therefore 

close to the asymptotic rate. Thus, in this application, the asymptotic 
rate of convergence is achieved when the mean errors are between one 
and five percent. If the initial guess is within (say) five percent 
of the solution, the performance of forty or fifty iterations will 
reduce the error by only one-half. The final values may have no more 
significant figures than the first guess. 

It is sometimes suggested that, in a numerical integration of the 
equations for fluid motion, if the time step is kept small, and if the 
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previous values of the stream function are used as a first guess for 
the new values, only a few iterations are needed for the solution of 
the Poisson equation. It is evident that this is not true if the 
number of mesh points is large. Even if only a very small change in 
the stream function, or pressure, is required, only a very small 
fraction of that change will be accomplished with a few iterations. 
Significant errors will be introduced at each time step. These errors, 
due to incomplete convergence, unlike round-off errors, are systematic 
in some sense and unacceptable in a long-term numerical integration. 

The errors in the direct solution by the Fourier transform arise 
from round-off and, in every case tested, were less by at least a 
factor of ten than the best iterative solution. The computing time 
was far less, and generally equal to that required for about nine 
iterations by the Over-relaxation. In the case of the 201 by 64 point 
array the Transform method produced a solution with errors less, by a 
factor of ten, than the best iterative solution, and reduced the com- 
puting time by a factor of twenty-five. 

E . DESCRIPTION OF THE PROGRAM AND COMPUTATIONAL PROCEDURE 

A listing of the computer program is reproduced as Appendix B. 

The names and functions of the various subroutines, the definitions 
of the parameters which control the operation of the program, as well 
as the dimensions of the various arrays, are given at the beginning of 
the calling program. The program is written in standard FORTRAN IV. 
Calculations were performed on an I.B.M. 360/67 computer at the W. R. 
Church computer center of the Naval Postgraduate School. Object codes 
were obtained with the I.B.M. H-level compiler, release 18. Data 
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storage requirements are the equivalent of 4 H N + 26 N + 13 M single- 
precision words, where N is the number of mesh points in the Y direction 

and M the number in the X direction. The. program itself requires the 

4 

equivalent of about 10 s ingle -prec is ion words of main storage. The 
calculations with a 201 by 64 point mesh required 13 seconds of com- 
puter time for the performance of a time step. The program incorporates 
provisions which enable a long integration to be performed in smaller 
segments. After a specified number of time steps all quantities in- 
volved in the integration are written out on a magnetic tape or disc. 

The integration can be restarted by reading the data from the tape or 
disc . 

The finite-difference mesh is specified by assigning, values to the 

parameters N, M and ALPHA. The last quantity is the basic period of 

the motion in units of L, one -ha If the distance between the parallel 

walls. It fixes the external streamwise dimension of the finite- 

difference mesh. In all calculations thus far ALPHA was equal to the 

wave length of the primary unstable mode. 

The time step, St, is controlled by the parameters FRACEL and 

XIAMDA. FRACEL is generally the operating parameter. Tt is equivalent 

to f in equation (4.43) and expresses the constraint on St associated 

with the advective instability of the central time differencing. 

XIAMDA is an upper limit on the ratio 

4 S t 
R 6 x 6 y* 

Fixing this ratio insures that the consistency condition, (4.30), is 
satisfied. All calculations were performed using a value of 1/4 for 
XIAMDA. 
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The parameter KINT specifies the number of time steps between the 
recalculations of 6t. At intervals of KINT time steps: the maximum 

components of horizontal and vertical velocity are determined; the 
value of 6t is adjusted; values of the vorticity are averaged over 
the last two steps; values of the stream function recalculated; and 
finally a single time step is performed by the Modified Euler method. 

Subroutines SETUP and START begin the operation. Values of the 
parameters associated with the finite-difference mesh and various 
other fixed quantities are established in the former. The initial 
values of the disturbance stream function and vorticity are establishes 
in the latter routine. 

The subroutines which perform the subsequent integration and their 
functions are listed below. 



ADVANC Time step by central differencing. 

STEP Time step by Modified Euler differencing. 

PRESUR Solution of Poisson's equation for the 

stream function. Calculation of 
vorticity at the walls. 

DIIARM Calculation of Discrete Fourier transform. 

TRISOL Solution of tridiagonal system of 

equations . 

TIMER Calculation of 6t . 



Subroutines LOAD and RELOAD perform the read and write operations 
associated with the interruption and restart of the integration. If 
the data parameter NPROB is less than zero then the main program calls 
subroutines SETUP and START for the initiation of a new problem. Upon 
the completion of the specified number of time steps the values of the 
various quantities are written on a tape or disc by subroutine RELOAD. 
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The remaining subroutines are associated with the calculation and 
printing of various derived quantities. The frequency at which the 
derived quantities are calculated and printed is controlled by the 
parameters 11SCIP and ISCIP. At intervals of MSCIP time steps the 
main program calls subroutine PRINT. Subroutines EDDYS , ENERGY, and 
STRESS are called in sequence. The current values of time, 6t, the 
normalized integrals of kinetic energy in the mean flow, the turbulent 
motion and of total kinetic energy are calculated and printed. At 
intervals of ISCIP time steps the main program calls subroutine PLOT. 
Values of the velocity and vorticity of the mea'n flow as well. *as the 
energy and Reynolds stress associated with the turbulent motion, as 
functions of Y, are printed and plotted. In addition three calling 
sequences are started in subroutine PIIESUR which involve subroutines 
SPCTRM, MODES and PLOTS P. Values of the spectral densities of squared 
vorticity and kinetic energy as functions of Y, as well as the normalized 
sums over Y of those quantities are calculated, printed and plotted. 

The functions of Y associated with specified disturbance modes are 
printed . 

F. FINITE DIFFERENCE REPRESENTATION OF THE LINEAR EIGENVALUE PROBLEM 
The stream function and vorticity are represented by functions 
which have the form of equation (3.12), but which are discrete in y. 
Derivatives with respect to y in equation (3.15) are replaced by the 
same second-order central differences employed for the non-linear case. 
The result is an algebraic eigenvalue problem. Equation (3.15) is 
approximated by equations of the form 

W? = i0$ 
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where, in the general case, W is an (N-2) x (N-2) complex matrix with 
eigenvalues 3 and eigenvectors $. Of particular interest is the effect 
of N, the number of mesh points between the walls, on the principle 
eigenvalues of W. 

Let 



= a + s* e i(ox-3*t) 



i(f = $ e 

J J 



r = 0 e -i(ax-pt) + 0>v e i(cxx-p*t) 
•J J 



(4.57) 



Equation (3.13) is satisfied if 

r 1 V, " ( 2 +a 2 6y 2 ) e T + e T , 

iaR [yp/a-u ) Gj + u"4>J - 2 — 

by 



0 (4 . 58) 



J = 2,3, • • • N- 1 



where from equation (3.14) 
$ 



6 J = 



J+l ~ (2 + a2&y2) 6 J + 9 J-1 



by' 



J = 2,3, • • ' N- 1 (4.59) 



The boundary conditions require that 



*i ' % = 0 



and that, in accordance with the previous formulation of the vorticity 
boundary conditions, 



where for 
equations 
given by 



g n = 



P 1 $ N-1 + P 2 $ N-2 + P 3 ^N-3 + P 4 ^N-4 



6y 



6 1 = 



P 1 ^2 + p 2 S + p 3 ^4 + p 4 *5 

by 2 



the various approximation for the boundary conditions, 



(4.12), (4.13), (4.14) and (4.15), p R , K = 1, 2, 3, 4, 



(4.60) 



are 
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TABLE III 



COEFFICIENTS FOR VORTICITY BOUNDARY CONDITIONS 





p i 


p 2 


P 3 


p 4 


first crcler 


2 


0 


0 


0 


second order 


4 


- 1/2 


0 


0 


third order 


6 


-3/2 


2/9 


0 


fourth order 


122/9 


-3 


. 8/9 


- 1/8 



Equations (4.58) and (4.59) can be written 



igGj = ia 



e T+l +CD0 T +e T1 

u e - u” $ + — - ^ — 

L J J JJ R6y 



9 J = 



$ + co$ + $ 

J+l J J-l 

6y 2 



where 

co = - (2 + 6y^) 

We define the following Matrices and vectors. 



(4.61) 



(4.62) 



(4.63) 



B = 


to 1 
1 co 1 
1 co 1 


P = 


Pi P 2 P 3 


p 4 




1 co 1 
1 co 

L 




p 4 


P 3 P 2 



(4.64) 



82 




(4.65) 



And adopt the notation for a diagonal matrix with diagonal elements 



Tl- 



J = 2 ,3 , • • -N-l 



A <V 



^3 

^4 



^N-2 

%-l 



(4.66) 



Equations (4.61) and (4.62) become 



if30 = A(iaU ) 0 - A (iaU ,f ) $ + 

j 




[B G + y) 



(4.67) 



0 = —■ B§ (4.68) 

6y 

The vorticity boundary conditions, equations (4.60), become 

7 = (4.69) 

6y 

The boundary conditions, = 0 are satisfied by equations (4.67) 

and (4.68). Substituting equations (4.68) and (4.69) in (4.67) we 
obtain 
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Multiplying through by B we obtain 




Thus 



ip$ - W$ 



(4.71) 



Where 




(4.72) 



For the linear equations the problem is simply the determination of 
the eigenvalues and eigenvectors of a complex matrix. Because of 
symmetry of the laminar velocity distribution the symmetric and anti- 
symmetric solutions of both the differential and algebraic equations 
are uncoupled. It is useful to take advantage of this fact to reduce 
the order of the matrix W. The well-known unstable solutions for the 
Plane Poiseuille flow are all symmetric. It appears that all the anti- 
symmetric solutions are stable (Grosch 1968). In the present under- 
taking only symmetric modes are considered. 

Thus we consider only the region 



and add the boundary conditions at the centerline (y = 0 , J = 



-1 < y < 0 



N+l 



• = 0 



In terms of the finite differences this is equivalent to 
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N+3 

2 



N-l _ N+3 

2 2 



$ 



N- 1 
2 



The matrix formulation remains the same except that the order of the 

N-l 

various matrices and vectors is reduced from N-2 to an< ^ equations 

(4.64) and (4.65) for B, P, §, 9 and y become 



’a) 1 


1 p l p 2 p 3 P 4 


1 to 1 


P - ! 


1 CO 1 


} 

t 

J 


1 to 1 


i 

! 

j 

j 


2 to 


j 





1 



(4.73) 




(4.74) 



The eigenvalues and eigenvectors of equation (4,71) were computed 
numerically by means of the Q R algorithm (Wilkinson 1965, Parlette 
1966). This algorithm, for the solution of an algebraic eigenvalue 
problem, was developed by Francis (1962). It involves a sequence of 
unitary transformations which yield an upper canonical form of the 
eigenvalue problem. The inverse of the matrix B in equation (4.70) 
is obtained by standard Gauss -Eliminat ion. All calculations were per- 
formed in double precision (sixteen decimal digits) on an I.B.M. 360/67 
computer. Standard I.B.M. scientific library subroutines, with some 
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minor modifications, were employed for the calculation of the inverse 
of B and the eigenvalues and eigenvectors of W. A listing of the 
FORTRAN program is reproduced as Appendix C. 
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V. RESULTS AND DISCUSSION 



A. LINEAR CALCULATIONS 

Symmetric solutions of the algebraic eigenvalue problem, equation 
(4.71), were obtained at regularly spaced points in a portion of the 
wave number -Reynolds number domain for various values of N ranging 
from 41 to 201. The investigation was confined to the vicinity of 
the unstable region at Reynolds numbers out to 25000. The calcu- 
lations yield, in each case, N/2 symmetric eigenfunctions and the N/2 
associated eigenvalues. A summary of selected results is depicted in 
Figs. 4-11. 

Preliminary calculations which employed the various order formu- 
lations of the vorticity boundary conditions, Table III, indicated 
substantial differences in the solutions with the first and second 
order formulations but only slight differences between solutions based 
on second and higher order formulations. Consequently the consistent 
second-order formulation was adopted for the subsequent linear calcu- 
lations, as well as for the numerical integration. 

Calculations with N = 41 revealed no indication of linear instability, 
(3^ < 0), at any point in the region investigated. Calculations were 
performed for wave numbers from .66 to .94 in increments of .02 for a 
regular sequence of Reynolds numbers out to 25000. In every case the 
calculated values of 3^ reached a local minimum at a Reynolds number 
less than 20000. Consequently it is believed that a second-order finite- 
difference representation of the equations of motion with N < 41, 

(6y > *05), completely obscures the phenomenon of linear instability. 

For N > 41, points on the neutral curve were obtained by interpolation. 
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Neutral curves, determined in this manner, for N = 61, 101 and 201, 
are shown in Fig. 4, along with the neutral curve established by 
Grosch. The agreement with established results is poor for N = 61, 
fair for N = 101 and for N = 201 reasonably good. Although it is not 
shown, the neutral curve for N = 161 differs only slightly from that 
for N = 201 and only on the upper branch. The accuracy of the repre- 
sentation, for a given value of N, is greatest in the vicinity of the 
lower branch of the neutral curve. This is illustrated in Fig. 5, 
where the calculated growth rates, 3^, are plotted versus the wave 
number, a, for a Reynolds number of 15000. Solutions for the eigen- 
values of the least stable mode, in a sequence of calculations with 
increasing N converge most rapidly for the smaller values of a. The 
same kind of behavior characterizes the convergence of the asymptotic 
expansions for the solution to the linear differential equation 
(Potter 1966) . 

The sets of calculated eigenvalues for the cases a - 1.0, 

R = 7333.3, 11333.3 and 25000 are plotted in Figs. 6, 7, and 8. The 
orderly arrangements of the eigenvalues in the complex plane are quite 
striking. None of the eigenvalues has a phase velocity, 3py greater 
than 1.5, the maximum velocity of the laminar flow, and, at the lower 
Reynolds numbers, the majority have phase speeds equal to the mean 
laminar velocity. The effect of increased Reynolds number is evidently 
the greatest on the most stable modes. As the Reynolds number is 
increased the eigenvalues approach the neutral line (3^. = 0). As 
R oo, the matrix W in equation (4.71) becomes Hermitian, In the limit 
all the solutions are neutrally stable. 
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A measure of the accuracy of the calculations can be obtained by 
comparison with selected results of Grosch (1968) and Thomas (1953). 
Calculated growth rates in the unstable region for N = 101, 1 6 1 and 
201, as well as values calculated by Grosch are plotted in Figs. 9 and 
10, for Reynolds numbers 7333.3 and 11333.3 respectively. The very 
slow convergence of the second order f inite-dif ference scheme in the 
vicinity of the upper branch of the neutral curve is evident. 

In terms of the reference units employed here , the primary eigen- 
values calculated by Thomas were, for a =■ 1.0, R = 1066.6 and 
R = 6666,7 respectively 

3 = .4847 + i .0393 
3 = .3563 - i .0055 

For these cases the present calculations with N = 201 yield 

3 = .4891 + i .0393 
3 = .3593 - i .0041 

A comparison of the calculated eigenvectors for the case R = 6666.7, 
a = 1.0 is given in Fig. 11. Results for the real part of the eigen- 
vector are identical to four significant digits. There are however 
small differences in the imaginary parts. Thomas employed a trans- 
formation which made his finite-difference approximation basically of 
fourth order. It is not surprising that a second-order scheme might 
produce slightly different results for a function with large derivatives 
and curvature. However test calculations, with a linearized form of 
the numerical integration, which employed the calculated linear » 
solutions as initial conditions, suggest that most of the differences 
are attributable to slight numerical errors in the present solutions. 
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FIGURE 5 EFFECT OF N OH CALCULATED GROWTH RATE 




CALCULATED EIGENVALUES - SYMMETRIC MODES 
R = 6666.7 a "1.0 N =20! 
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CALCULATED EIGENVALUES - SYMMETRIC MODES 
R - 1 13334 ^ = 1.0 N = 20! 
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CALCULATED EIGENVALUES- SYMMETRIC MODES 
R = 25,000 * = 1.0 N = 201 
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ISURE 10 COMPARISON OF CALCULATED GROWTH RATES. 



0.020 
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The errors associated with a particular eigen-solution by the Q. R. 
algorithm are dependent, among other things, on the ratio of the 
modulus of the associated eigenvalue to the largest such modulus. 
Decreased values of that ratio imply increased errors in the solution. 
In the present application the eigenvalue associated witli the least 
stable mode has the smallest modulus, and it is reasonable to expect 
that there might be some accumulation of numerical error in the 
solution for that mode. 

In retrospect it seems that, for the linear solutions, an iterative' 
method which focuses on a single mode, such as the approach introduced 
by Nachtscheim (1964), might have been more appropriate for the limited 
and definite purposes of the present undertaking. However in any 
application in which there is interest in the higher modes the matrix 
approach has definite advantages. 

B . NUMERICAL INTEGRATIONS 

I . Pre liminary 

The linear calculations, in addition to providing a measure of 
the resolution required to represent adequately the two dimensional 
motions, provide a rational means of establishing initial conditions 
for the numerical integration. In particular they provide a unique 
means of testing the stability (in the linear range) and accuracy of 
the finite-difference formulation and solution. The calculation of 
the non-linear terms in the finite-difference solution can be bypassed 
by minor changes in the FORTRAN code. If this is done, the scheme 
should be capable of reproducing the calculated linear eigensolutions . 
Specifically, if an unstable mode is employed to establish the initial 
conditions, the numerical integration should produce the predicted 
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growth rate and phase speed. The variation of the calculated turbulent: 

-20 t 

kinetic energy with time should be, from equation (3.18), e H x . 
Calculations of the spectral density of turbulent kinetic energy and 
-squared vorticity should reveal no spurious excitation in other 
represented wave lengths. Finally, perhaps the most stringent require- 
ment is that the performance of a large number of time steps should noc 
produce any significant distortion of the. normalized eigenvector. 

If the scheme can be demonstrated to be stable and accurate, at 
least in the linear features, then, since the virtues of the particular 
non-linear representation have been fairly well established elsewhere, 
calculations can be performed in the non-linear range with confidence 
in tiie accuracy of the results. 

The resolution requirements indicated by the linear calculations 
are very great: about 201 mesh points are required in the Y direction. 
The requirements for resolution in the X direction are determined by 
the extent to which higher harmonics, or smaller scales of motion, are 
excited by the non-linear terms. M points in the X direction represent 
M/2 Fourier modes. As a result of the well-known constraints on 
energy cascades in a two-dimensional motion the extent of the excitation 
of the higher harmonics should not be great. Kraichnan (1967) has pre- 
dicted a -3 power law for the spectral density of turbulent kinetic 
energy in the smaller scales. If this is the case, then the energy in 
wave number 10, for example, should be less by three decades than that 
in the primary mode of the disturbance. If the basic period of the 
motion, the external downstream dimension of the finite -d if ference mesh, 
is chosen to be equal to the period of the unstable mode, a mesh with 
M = 64 , which represents 32 Fourier modes, should be adequate. 
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2. Results 



A tabulation of information which serves to identify the 
various integrations or runs which will be referred to in this section 
is given in Table IV. 

TABLE IV 

SUMMARY OF NUMERICAL INTEGRATIONS 



Run 


Initial 
Con- 
dit ions 


2a 


E(o) 


Number 

of 

Time Steps 


Time 


Computer 

Time 

(Iirs) 


i 


B* 


.3 


.4xl0 -1 


1000 

1 


11.7 


00 

L_ 


2 


A* 


.L/5 


1 .25x10"'' 


1550 


20.0 


6.1 


3 


A 


.01 


.5xl0" 4 


800 


8.3 


3.8 


4 


A 


. 01/10 


.5xl0' 3 ' 800 


8.3 


3.8 


5 




. 1 


I . 5xl0 -2 


800 


8.2 


3.8 


6 


A 


. l./Io 


. 5x10” 1 


800 


7.6 


' j 

3.8 


7 


A 


1.0 


.5 


3700 


11.0 


17.5 


8 


A 


. 1/5 


. 25x 10” 1 


5900 


61.2 


28.0 


9 


A 


. Ia/2 


. 25xl0" 2 


2700 


52.5 


12.0 j 



Initial conditions: 

■ o f, ioix , -iax 0 

T “ a / $ T e K + $* e K > 

R , J C J J ^ 

A - R = 6666.7, a = 1.0 

§ from Thomas (1953) 
j 

B - R = 15000, a = .86 

$ from present calculations 
J 

*non- linear terms suppressed. 
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As indicated in the table, runs one and two were performed 
with the non-linear terms suppressed. The remaining runs included 
the non-linear interactions and were at a Reynolds number of 6666.7. 

The linear solution calculated by Thomas was employed as an initial 
condition. They differ in the initial amplitude, a. 

Preliminary test calculations, with the non-linear terms 
suppressed, which employed, as initial conditions, the primary eigen- 
functions calculated here indicated not only the accuracy of the 
integration scheme, but also a slight inaccuracy in the calculated 
eigenvectors. Run number one is typical of these calculations, most 
of which were performed in fairly coarse arrays. The logarithm of the 
ratio of the turbulent energy to its initial value, in that run, is 
plotted as a function of time, along with predicted rate in Fig. 12. 
During the course of an initial transient which persists for about six 
units of time, the growth rate increases and then decreases to about 
the predicted value. Although the errors in the growth rate, during 
that period, are large the asymptotic value is very good. Other calcu-^ 
lations with N = 41 and a decaying disturbance exhibited the same 
behavior. In both cases the predicted rate was achieved and persisted 
after the transient. Reference to equation (3.19) suggests the cause 
of the transient, the presence of small errors in the calculated eigen- 
vectors. In effect, small contributions from higher modes are present 
in the linear solution for the least stable mode. Figures 6 and 7 
illustrate the fact that the stability of the great majority of the 
higher linear modes far exceeds the instability of the single unstable 
mode; there are wide variations in the exponential rates of growth or 
decay. After several units of time the energy associated with the 
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errors is washed away. The integration routine, if it is accurate, 
acts as a filter and tends to suppress contributions from all but the 
least stable mode when the non-linear terms are omitted. 

A subsequent run, (run number two), at a Reynolds number of 
6666.7 which employed the eigenvector calculated by Thomas as initial 
conditions tends to confirm this suggestion. For this run, the 
logarithm of the energy ratio, along with both the growth rates pre- 
dicted in the present calculations and that of Thomas are plotted in 
Fig. 13. There is a very slow decrease in the growth rate out to about 
ten units of time. After that point the variation is linear and very 
close to the growth rate predicted by the second-order calculations. 

In effect, the eigenvector calculated by Thomas in a basically fourth- 
order formulation is a better approximation of the second order 
solution than those obtained in the present work. After 1000 time 
steps the vector was renormalized and printed, in order to obtain a 
measure of the distortion. The variations, where they occur, are 
plotted in Fig. 14. The maximum distortion in amplitude was .035/' 
and in phase 1.5$. 

Calculated spectral densities of kinetic energy and squared 
vorticity revealed no indication of spurious excitation of other 
spatial wave lengths in any of the linear runs. 

Thus it may be concluded that the integration scheme is stable 
and highly accurate. 

Although the performance of the integration routine in the 
linear tests was most satisfying, the inaccuracies in the calculated 
eigenvectors and the resulting transient in the linear integration 
inhibit a most useful application of the program. That is, to the 
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determination of critical amplitudes at sub-critical values of the 
Reynolds number. For this application the representation of the 
linear behavior in the early stages should be more accurate than that 
which results from the present calculation. However, more accurate 
methods, for the determination of the solutions of the algebraic 
equations in a consistent second-order model of the linear problem, 
are available. In particular the refinements of the iterative approach 
of Nachtsheim (1964), which were developed and employed by Reynolds are 
most promising. 

In the subsequent non-linear calculations at R = 6666.7 and 
a = 1.0, the linear solution, calculated by Thomas (1953), was employed 
for the initial condition. The amplitude factor, a, is based on the 
normalization of the linear solution given by Thomas. The objectives 
of the non-linear calculations were, to determine the effect of various 
amplitudes on the initial development of the disturbance, and, if 
possible, to achieve an equilibrium state for the two-dimensional 
turbulence . 

In runs three through six the initial amplitudes were such that 

the initial kinetic energy of the disturbance increased, by decades, 

-4 - 1 

from .5 x 10 to .5 x 10 . For reference, the normalized sum of 

kinetic energy in the laminar motion is . 6 in the present units. The 
growth of turbulent energy for these calculations, as well as for run 
number nine, are shown in Fig. 15. For run number three, no significant 
departures from the linear behavior were observed until after five units 
of time. This provides some measure of the amplitude at which the non- 
linear interactions become significant. The limit in the growth of the 
turbulent energy in run number six suggests the existence of an 
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equilibrium state, in terms of kinetic energy, in the range 

.005 < E < .05. 

In order to verify the limiting behavior of the turbulence a 
long calculation was performed at an even higher energy level. This 
was run number five. In this case the initial kinetic energy of the 
disturbance was of the same order as the energy in the laminar flow. 

The subsequent variations in the kinetic energies of the mean flow, 
the turbulent motion, and of the energy in the primary mode of the 
disturbance, as well as the total kinetic energy in the domain are 
shown in Fig. 16. All quantities are normalized with respect to the 
initial value of the total energy. Points are plotted at intervals 
of 100 time steps. It is the variation of the total energy which most 
clearly indicates the trend. It decayed monotonica lly , and by the end 
of the run about 13/a of the kinetic energy in the system had been 
dissipated. The early behavior of the turbulence is dominated by 
large exchanges of energy between it and the mean flow. 

Run number eight, in which the initial energy of the disturbance 
was .025, produced an approximate equilibrium condition. The variation 
of the various energy terms is shown in Fig. 17. After a slight decrease, 
the total energy remained virtually constant for the last 25 units of 
time. Again, there is a periodic exchange of energy between the mean 
flow and the disturbance. However in this case the period of the 
oscillation was more than five times as great as in the previous case, 
and, in the course of the development, the amplitude decreased markedly. 

A comparison of the variations of the total kinetic energy in the two- 
dimensional turbulence to the energy in the primary mode suggests that 
the oscillations are associated with transients in the development of 
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the higher wave numbers. After 47 units of time the amplitude of the 
oscillation was much reduced , and the difference between the two 
kinetic energy terms, which is the total kinetic energy in the higher 
wave numbers, remained nearly constant. Mean velocity profiles at 
extreme points in the early variations of the disturbances are given 
in Fig. 19. Plots of the mean velocity, mean vorticity and turbulent 
kinetic energy near the final state are given in Fig. 20. The mean 
velocity profile is only slightly distorted although near the walls 
there are significant changes in its curvature. Calculated spectral 
densities of the turbulent kinetic energy are consistent with the -3 
power law proposed by Kraichnan. Values of the spectral density of 
kinetic energy at y = + . 85 are given in Fig. 21, mean values over y, 
in Fig . 22 . 

The function of y which describes the primary mode remains 
symmetric. In addition, it turns out that modes three and five are 
also symmetric while modes two and four are anti-symmetric. Although 
the remaining mode shapes were not determined in the calculations, it 
appears likely that the odd-numbered modes are symmetric and the even- 
numbered modes ant i -symmetric . Mode shapes, in terms of amplitude and 
phase, for modes one through four are given in Figs. 24 through 27. 

The initial shape for mode one is given in Fig. 23. It is interesting 
to note that, although the symmetry and anti-symmetry of the odd and 
even-numbered modes does not necessarily follow from equations (3.4), 
it does follow from the equations treated by Watson (I960) and Reynolds 
(1967), in which certain of the non-linear terms are ignored. 

The last calculation, run number nine, was started with a 
disturbance energy, less by a factor of ten, than that of run eight. 
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In this case the variations of the various energy terms were nearly 
monotonic. After 50 units of time the disturbance energy had increased 
by more than a factor of three. It appears likely that if this calcu- 
lation were continued the equilibrium energy level would have been 
attained without the large oscillations which characterized the initia] 
behavior of the previous run. 



106 



I 



T 1 r 



T 



T — — *1 — 




107 



FIGURE 12 TEST FOR UNEAR SOLUTION . 
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SGURE 14 RUN NO. 2 DISTORTION OF EIGENVECTOR AFTER 1000 STEPS. 
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Tims Dimensionless Units 
FIGURE 15 

EFFECT OF AMPLITUDE ON INITIAL GROWTH RATE 




FIGURE IS ENERGIES VS. TIME - RUN 7 



in 




FIGURE I7A ENERGIES VS. TIME- RUN 8 
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FIGURE I7C ENERGIES VS. TIME - RUNS Cont. 
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.FIGURE 13 MEAN VELOCITY PROFILES , RUN 8 
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FIGURE 20. RUN 8. TIME =53.2 
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FIGURE 21 
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FIGURE 22 
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VI. CONCLUSIONS 



The f inite -d if ference representation of the vorticity transport 
equations, and the explicit treatment of the boundary values of the 
vorticity are consistent, stable and accurate. This has been demon- 
strated, in the linear range, by comparison with the established 
results of the linear theory. 

The application of the discrete Fourier transform to the solution 
of the system of equations, associated with the representation of the 
Poisson equation, has produced significant improvements in both the 
speed and the accuracy of the calculations. Numerical tests indicate 
that the usual iterative schemes are inadequate if the number of mesh 
points is large. 

Extensive calculations based on a consistent second-order repre- 
sentation of the linear eigenvalue problem provided a measure of the 
effect of mesh size on the accuracy of the representation of the large- 
scale unstable motions. The results indicate that the mesh sizes 
employed in previous investigations were often inadequate to represent 
linear instabilities . 

Calculations in the non-linear range, at a moderate Reynolds number, 
indicate the existence of an equilibrium amplitude of the disturbance 
motion. In a long-term numerical integration, an approximate value of 
the equilibrium amplitude has been obtained. At the equilibrium state, 
the mean intensity of the two-dimensional turbulence is somewhat greater 
than measured values in real turbulent flows. This reflects the reduced 
dissipation due to the constraints on energy cascades in two-dimensional 
turbulence. The calculated spectral densities of kinetic energy are 
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consistent with the -3 power law predicted by Kraichnan for the two- 
dimensional turbulence. 

The present program provides a means, of demonstrated accuracy, 
for the extensive investigation of finite -amplitude instability and 
two-dimensional turbulence . 
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VII. RECOMMENDATIONS FOR FUTURE WORK 



The computer programs developed in the course of this investigation 
for both the solution of the linear eigenvalue problem and the numerical 
integration have yet to be fully exploited. A number of interesting 
applications remain. Among them are: 

1. More.' extensive calculations for the solutions of the linear 
eigenvalue problem, in particular for the asymmetric modes. 

This will require only minor modifications of the program. 

2. The development of an iterative procedure for the refinement 
of the solution for the eigenvector associated with the least 
stable mode. This will facilitate the investigation of finite - 
amplitude instabilities at sub-critical values of the Reynolds 
number . 

3. A continuation of the integration for run number nine. The 
results to date show the disturbance increasing smoothly and 
monotonica lly at a decreasing rate. If, as seems likely, the 
energy levels off at the values associated with run eight, the 
indication of the equilibrium state will be more emphatic. 

4. A more extensive invest igat ion of the structure of the two- 
dimensional turbulent motion. This could be useful for the 
testing and development: of heuristic models, in particular, 
the promising model developed by Gawain and Pritchett (1970). 

5. Finally, some efforts should be made to further improve the 
speed of the calculations. A computer analysis of the operations 
in the program for the numerical integration indicates that sub- 
stantial improvements in speed could be achieved by coding just 
Subroutine ADVANC in assembly language. 
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APPENDIX A 



DERIVATION OF DISCRETE RELATIONS FOR THE SPECTRAL DENSITIES 
OF SQUARED VORTICITY AND KINETIC ENERGY 

Associated with discrete Fourier transform pair, equations (4.46), 

is the convolution 



M . M 

1 2 T T* e' 1 °n,s = S F F* .. (A.l) 

M - n n Tr , I< K-s+1 

n=l K=j- 

where 0 is as given by equation (4.47). Equation (A.l) can be 
n > s 

derived by substituting the appropriate transform or inverse relation, 
and reversing the order of the summation. 

For s = 0 ve have the spectral density relation 



1 



M 



M 



- ^ T T* - X F v F* 

M - n n - K K 
n=l K=1 



(A. 2) 



The mean value of squared vorticity is given by 

a j = M ^ r K,J r K, J 
K— i 



(A. 3) 



Then from equations (A. 3), (4.49) and (A. 2), 



i M 

a . = — 2 G G* _ 

j .2 . n , J n , J 

J M n=l 



Since F is real, 
K , J 



1 C 2 M ^ 2 

= -+r 2 g: _ + 2 2 G 



G* 



“J ,.2 £ 1,J n _2 ~n,J ”n,J 



M 



Si/ 2+1, J 1 

W/ 



(A. 4) 



Thus the quantities G. , 2G G* , 2G„ , G* etc. are the 

T j J ^ j vj Z j J J ^ J J ^ J 

spectral components of mean squared vorticity. Notice that, in the 
present notation the first spectral term, n=l, refers to the contribution 
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of the disturbance to the mean flow. If the laminar flow vorticity 
is included then the first term becomes 



Of interest also are 
a^ over y, 



< g i,j + 3 ^>‘ 



the spectral components 



of the normalized sum of 



1 

N- 1 




N-l 

E 

J=2 



a j + 



} 



which are obtained by performing the same operation on the various 
terms on the right of equation (A. 3). 

Expressions for the spectral components of kinetic energy, at a 
given value of J, or their normalized sums over y, can be derived in 
a similar way. The components of turbulent velocity can be expressed 
as 



K, J 



1 M 

r. * U e n,K 
M - n , J 
n=l 



1 



M 



(A. 5) 



K,J 



77 S V T 
M - n , J 
n~ i 



e 1 9 n,K 



where, from equations (2.25) and (4.49) U , and V are given by 

n , J n , j 



n, J 



F - F 

n,J+l n , J - 1 

26y 

i s in 9 



(A. 6) 



V 



n, J 



6x 



f 



n , J 



Thus from equations (4.24) and (A. 2) the mean value of turbulent 
kinetic energy is given by 
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A 




2M 



M 

E 

n=2 



U 



n , J 




+ 



1_ 

2M 



M 

E 

n-2 





Since u and v are real 
1C } 'J y J 



A i f A x 2 2 l 

E=-z~/}2 E ( U _ U* _ + V — V* ) + U , + V , ( 

J 2M C n =2 ^ n,J n>J n,J n ’ Jy ^ M/2+1, J M/2+1, JJ 



(A. 7) 



Finally, the spectral components for 



A 

E 



N-l 



N-l 

E 

J=2 



-X 




are obtained by performing the same operation on the terms on the 
right in equation (A. 7). 
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APPENDIX B 



LISTING OF THE FORTRAN CODE FOR THE NUMERICAL INTEGRATION 
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APPENDIX C 



LISTING OF THE FORTRAN CODE FOR THE SOLUTION 
OF THE ALGEBRAIC EIGENVALUE PROBLEM 
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